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^A  normal-mode  program  for  a sound-speed  profile  of  an  arbitrary  number  of  layers  has  been  constructed. 

It  has  been  used  extensively  and  successfully  for  12  years  to  compute  sound  propagation  in  idealized  underwater 
acoustic  ducts.  Documentation  is  contained  for  those  who  wish  to  use  or  modify  this  program.  The  FORTRAN 
statements  are  given  for  both  the  normal-mode  and  a mode-follower  program.  The  numerical  analysis  necessary 
for  computing  modified  Hankel  functions  of  order  1/3  is  given.  The  analysis  includes  a continued  fraction 
technique,  x 
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OBJECTIVE 

Translate  well-known  differential  equation  solutions  into  a working  program  to  com- 
pute propagation  in  underwater  acoustic  ducts.  Document  the  program  methods,  to  assist 
users  of  this  and  similar  programs. 


RESULTS 

1 . An  effective  program  for  computing  propagation  loss  in  a layered  ocean  by  normal 
modes  has  been  developed.  Complete  documentation  for  the  program  is  contained  herein. 

2.  Sediment  layers  are  modeled  as  fluids  in  which  densities,  sound  speeds,  and  absorp- 
tion can  be  specified.  This  permits  a complete  wave  solution  for  bottom  reflected  sound 
energy. 

3.  A continued  fraction  technique  for  evaluating  asymptotic  series  is  shown  to  give 
superior  results  in  evaluating  the  auxiliary  functions  required  in  this  program,  the  modified 
Hankel  functions  of  order  1/3. 

4.  A mode  follower  program  given  here  is  useful  in  tracing  eigenvalues.  Such  traces 
are  needed  to  understand  the  eigenvalue  structure. 

RECOMMENDATIONS 

1.  Improve  the  mode  locating  ability  of  this  normal-mode  program  to  make  it  self- 
contained.  It  currently  requires  user  interaction  to  locate  eigenvalues. 

2.  Investigate  methods  to  incorporate  the  effect  of  rough  boundaries  into  this 
program. 
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INTRODUCTION 


This  report  describes  a normal-mode  program  that  has  been  used  successfully  for  12 
years  to  compute  sound  propagation  in  idealized  underwater  acoustic  ducts.  The  theory  and 
considerations  used  in  developing  the  program  are  discussed  here,  and  a copy  of  the  FOR- 
TRAN statements  are  included  as  appendix  A.  Appendix  B consists  of  sample  inputs  and 
outputs  to  assist  users  in  gaining  familiarity  with  the  program.  It  is  hoped  that  this  report 
contains  sufficient  information  to  allow  a user  to  run  the  program  and  to  modify  it  as  desired. 

This  program  follows  the  methods  developed  by  Furry  and  Freehoffer  (ref  I)  to  com- 
pute electromagnetic  propagation  in  the  1 940s.  Marsh  adapted  these  methods  to  underwater 
sound  in  his  doctoral  thesis  (ref  2).  Using  this  material,  Pedersen,  at  NOSC  in  the  late  1950s, 
adapted  the  method  to  digital  computers  and  developed  the  programs  to  compute  the  auxil- 
iary functions.  This  original  program  used  two  layers  to  define  the  sound-speed  profile  (ref  3). 
This  program  was  expanded  to  three  layers  by  DF  Gordon  and  RF  Hosmer  and  finally  to  the 
multiple-layer  program  reported  here.  In  this  program  the  only  constraints  on  the  number  of 
layers  are  computer  space  and  running  time.  The  program  is  normally  configured  to  permit 
up  to  1 2 layers. 

The  earlier  programs  were  used  to  study  sound  propagation  in  ocean  surface  ducts. 
Programs  that  permit  more  layers  have  proven  useful  also  for  studying  propagation  in  the 
deep  ocean,  although  the  number  of  modes  required  generally  limits  computations  to  fre- 
quencies below  300  Hz.  The  multiple-layer  program  has  also  proven  useful  in  modeling  sedi- 
ment layers  and  thus  in  computing  shallow-water  propagation. 

The  principal  limitation  in  the  application  of  this  program  to  real-world  situations  is 
the  requirement  of  ideal  conditions:  boundaries  must  be  smooth  and  horizontal,  and  no 
variation  of  boundary  conditions  with  range  is  permitted.  Despite  this  limitation,  the  pro- 
gram has  proven  useful  in  predicting  and  explaining  acoustic  propagation  and  has  applications 
in  a number  of  related  areas.  1'hese  include  checking  other  types  of  wave-theory  models  or 
corrections  such  as  caustic  corrections;  determining  group  velocities,  dispersion  curves,  and 
reflection  coefficients;  and  determining  acoustic  coupling  between  ducts. 

The  following  paragraphs  describe  the  specific  topics  covered  by  the  sections  in  this 
report.  In  GENERAL  SOLUTION  are  the  equations  required  to  solve  the  wave  equation 
with  the  boundary  conditions  used  here.  DETERMINANT  is  part  of  the  basic  solution  but 
is  concerned  with  the  particular  numerical  method  used  in  this  program  to  evaluate  the  con- 
ditions imposed  by  the  boundaries.  Other  approaches  could  be  used  instead.  A later  section, 
NUMERICAL  BREAKDOWN,  is  also  part  of  the  basic  solution,  but  deals  with  special  numer- 
ical problems  that  have  arisen  but  are  not  apparent  from  the  basic  equations. 


1 . The  Bilinear  Modlfled-Index  Profile,  by  WH  Furry,  in  Propagation  of  Short  Radio  Waves,  DE  Kerr,  ed; 
MIT  Rad  Lab  series,  vol  13,  p 140-168,  McGraw-Hill.  New  York,  1951. 

2.  Navy  Underwater  Sound  Laboratory  Report  111,  Theory  of  the  Anomalous  Propagation  of  Acoustic 
Waves  in  the  Ocean,  by  HW  Marsh,  1950. 

3.  Normal-Mode  Theory  Applied  to  Short-Range  Propagation  in  an  Underwater  Acoustic  Surface  Duct,  by 
MA  Pedersen  and  DF  Gordon;  J Acoust  Soc  Am,  vol  37,  p 105-1 1 8,  January  1965. 
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FINDING  EIGENVALUES  deals  with  the  philosophy  of  eigenvalue  location  employed 
by  this  program,  which  essentially  leaves  this  function  to  the  user,  the  program  only  serving  as 
a tool.  It  shows  how  the  program  is  used  to  make  computations. 

Several  “automatic”  mode  Finding  versions  of  this  program  have  been  developed  to 
the  point  of  accommodating  certain  classes  of  profiles.  However,  they  need  further  develop- 
ment and  have  not  yet  been  reported. 

SOUND  SPEED  PROFILE  indicates  the  required  equations  for  curve  fitting  and  the 
various  ways  the  sound  speed  can  be  read  in  on  cards.  A continuous  water  profile  can  be 
entered  quite  simply,  but  sediment  layers  with  sound  speed  discontinuities  and  absorption 
gradients  can  become  complicated. 

REFLECTION  COEFFICIENTS  AND  OTHER  AUXILIARY  OUTPUTS  describes  a 
short  subroutine  that  computes  reflection  coefficients  for  any  mode  at  a given  profile  inter- 
face. Intermode  interference  lengths  and  mode  damping  coefficients  are  also  discussed. 

COMPUTATION  OF  THE  MODIFIED  HANKEL  FUNCTIONS  gives  the  analysis 
necessary  for  computing  these  functions.  The  use  of  continued  fractions  to  evaluate  an 
asymptotic  series  is  discussed.  To  facilitate  running  the  program  on  computers  of  different 
word  length,  this  section  provides  the  information  required  to  optimize  the  functions  for 
the  different  word  lengths. 

MODE  FOLLOWER  PROGRAM  describes  a separate  but  related  program  for  investigating 
the  eigenvalues  themselves  rather  than  using  them  to  compute  propagation  losses. 


GENERAL  SOLUTION 

The  derivation  of  the  normal-mode  solution  has  been  discussed  from  various  points 
of  view  (eg  ref  1, 4,  5).  Only  an  outline  is  given  here.  In  general,  the  time-independent  wave 
equation  is  written  in  polar  coordinates  and  the  azimuthal  coordinate  is  dropped  under  the 
assumption  that  ths  field  is  independent  of  azimuthal  direction.  Thus 

(1/r)  (d/dr)  [r(d*/dr)l  + (d2*/3z2)  + (w2/c2)  * = 0,  (1) 

where  \p  is  the  velocity  potential,  c the  sound  speed,  and  the  independent  variables  are  depth, 
z,  and  range,  r. 

Equation  (1)  is  then  separated  into  range-  and  depth-dependent  parts  with  a separation 
constant  X.  The  separation  is  possible  when  the  sound  speed  is  a function  of  depth  only. 

After  accounting  for  the  source  discontinuity  and  the  outgoing  radiation  condition,  integrating 
over  all  real  values  of  the  separation  constant,  and  normalizing,  one  can  find  the  solution  for 
a field  point  in  terms  of  propagation  loss  H as  follows: 


4.  Naval  Air  Development  Center  Report  NADC-72002-AE,  Normal  Mode  Solutions  and  Computer  Programs 
for  Underwater  Sound  Propagation,  by  CL  Bartberger  and  LL  Ackler,  4 April  1973. 

5.  A Normal  Mode  Theory  of  an  Underwater  Acoustic  Duct  by  Means  of  Green’s  Functions,  by  RL 
Deavenport;  Radio  Sci,  vol  I,  p 709-724,  1966. 


+ aAr> 


(2) 


H* -10  log 


N 


psphir  y H‘(Xnr)Un(z)Un(z0) 
n=1 


2 

where  r is  the  range,  zq  is  the  source  depth,  z is  the  receiver  depth,  Hq  is  the  Hankel  function 
of  order  zero,  second  type,  Xn  is  the  nth  eigenvalue,  Un  is  the  depth  function  for  mode  n. 
and  ps  and  % are  the  densities  at  source  and  receiver.  The  sum  is  over  the  number  of  modes, 
N.  making  a significant  contribution.  The  final  term  contains  the  volume  attenuation  coeffi- 
cient, a A-  From  Thorp  (ref  6),  aA  in  dB/m  is  computed  by  the  relationship 


0.9144oA  = 0.0001  F2/(l  + F2)  + 0.04  F2/(4100  + F2), 


(31 


where  F is  the  frequency  in  kHz.  Improved  equations  or  those  for  specific  ocean  areas  can 
be  easily  substituted.  The  depth  function,  Un,  is  a solution  to  the  depth-dependent  part  of 
the  separated  wave  equation 

d2U/dz2  + (w2/c2(z)  - X2]  U = 0,  (4) 

where 


= 2>rf 

and  f is  the  frequency,  in  Hz. 

A closed-form  solution  to  eq  (4)  can  be  obtained  when  the  reciprocal  sound  speed 
squared  or  squared  index  of  refraction  is  a linear  function  of  depth.  That  form  is  used  in 
this  program,  and  sound  speed  in  each  layer  is  expressed  as  follows: 

|Cj/c(z)|2=l  • 2 7j  (z  - Zj)/cj,  (5) 

where  Cj,  zj,  and  7j  are  the  sound  speed,  depth,  and  sound-speed  gradient,  respectively,  at  the 
top  of  layer  i.  Up  to  1 2 such  layers  are  permitted  by  the  program,  for  modeling  the  sound- 
speed  profile. 

With  this  expression  for  sound  speed,  solutions  to  eq  (4)  can  be  expressed  in  terms 
of  solutions  to  Stoke*’  equation 

h”  + zh  = 0.  (6) 

Only  a simple  change  in  independent  variable  is  required  from  z to  f , where 

fj(z)  “ [a3(z  - zj)  + w2/c?  - X2]  /af  (7) 

and 

* -27j  co2/ct3.  (8) 


The  solutions  to  Stokes’  equation  that  are  used  are  the  modified  Hankel  functions  of  order 
1/3,  h j(f)  and  h2(f)-  The  depth  function  is  a linear  combination  of  these  two  independent 
solutions: 


Fn.i<I>*'Vjhl«n)tBr.4h2«n>-  <9> 

where  Fn  is  the  unnormalized  form  of  Un.  The  coefficients  an>j  and  Bn>j  for  mode  n in  layer 
i are  determined  to  satisfy  boundary  conditions,  which  will  be  listed  below.  Values  of  An  for 
which  the  boundary  conditions  can  be  satisfied  are  the  eigenvalues. 

The  first  boundary  condition  is  the  radiation  condition.  It  is  satisfied  by  using  a nega- 
tive sound-speed  gradient  in  the  deepest  layer,  which  extends  to  infinite  depth,  and  by  letting 
the  depth  function  there  be  proportional  to  h2  only.  That  is, 

Fn(z)  = Bnh2(rn).  (10) 

At  the  surface  the  depth  function  is  zero: 

F„(0) = 0,  (11) 

and  al  layer  interfaces,  pU  and  its  depth  derivative  are  continuous: 

PiFn,i<z)  = Pi+lFn,i+l<z>;  (12) 

dFni(z)/dz  = dFni+,(z)/dz.  (13) 

Here  pj  is  the  density  in  layer  i,  and  the  excess  acoustic  pressure,  p,  is  given  by 
p = pU. 

If  U is  assumed  to  be  the  vertical  component  of  the  velocity  potential,  eq  (12)  and  (13)  are 
equivalent  to  requiring  that  the  pressure  and  the  vertical  component  of  particle  velocity  be 
continuous  across  the  layer  interface. 


Applying  these  boundary  conditions  to  a sound-speed  profile  consisting  of  M layers 
results  in  2M  - 1 linear  equations  in  hj  and  f»2.  They  are  homogeneous  in  that  the  constant 
is  zero  in  each  equation.  There  are  M-l  coefficients  Aj  to  be  determined  and  M coefficients 
Bj.  These  coefficients  can  therefore  be  determined  within  a constant  of  proportionality  D, 
provided  the  system  of  equations  is  linearly  dependent.  That  is,  the  2M  - 1 square  matrix  of 
coefficients  of  Aj  and  Bj  must  be  of  rank  2M  - 2 or  less.  Its  determinant  will  then  be  zero. 
This  is  the  eigenvalue  condition.  Values  of  X must  be  found  which  make  the  determinant 
zero.  This  determinant,  G,  is  discussed  in  more  detail  in  a later  section. 

Zeroes  of  the  determinant,  G,  are  found  by  using  the  secant  method.  The  variable  in 
this  iterative  method  can  as  well  be  some  function  of  X as  X itself,  and  we  use  the  following 
complex  phase  velocity  (v): 

Xn  = w/vn.  04) 

To  find  a v that  is  a root  of  G requires  an  initial  guess,  v j , where  the  subscript  1 refers  to  the 
step  in  the  iteration  and  a small  increment,  8 j . Each  succeeding  estimate  is  given  by  the 
relationship 


VI  *vi  + 6r 


where 


6j--(vj-vj_,)Gj/(Gj-Gj_,>. 


The  details  of  this  iterative  process  are  given  in  a later  section. 

When  an  eigenvalue  vn  is  found,  the  coefficients  are  then  evaluated.  One  coefficient 
can  be  given  an  arbitrary  value,  so  A|  is  set  to  |(0)| . From  eq  (1 1).  B|  is  then 

-p I h | If  i(0)l . Pairs  of  equations  (eg  ( 12)  and  (13))  for  each  successive  interface  can  then 
be  used  to  evaluate  the  next  A;  and  Bj  as  discussed  later. 

Finally  the  normalizing  factor,  Dn,  for  mode  n is  obtained  by  the  relationship 


D„-  / 


This  equation  follows  from  the  orthogonality  of  the  depth  functions.  It  is  not  the  pressure, 
however,  which  is  proportional  to  pU,  but  p^U  that  is  orthogonal  (ref  7).  Therefore,  Dn 


however,  which  is  proportional  to  pU,  but  p'HJ  that  is 
must  be  determined  such  that  the  integral  of  pU-  is  1. 


From  Stokes’  equation  (eq  (6))  and  eq  (7-9),  the  integral  of  F“  takes  the  form 


zi+l  f / / ,"|  ZH 

J F*(f)dz  « fj(z)  F“(?)/tfj  + F'2($)h. 

zi  L ■*  zi 


Therefore 


II-  I 

Dn  * ~pfw2/ai  +£  Pi ( f 4(Zi+ 1 )/a|  - pjf i+ , (zi+  j )/(ai+ |Pi+1)l  F“(zi+1) 
i-1 

+ (pi/af-Pi+,/af+^F;2(zi+i)  , (18) 

where  eq  (1 2)  and  (13)  have  been  used  to  combine  terms  at  each  interface.  The  derivative  of 
F takes  the  form 

Fi(Vl)-«i(Aihiiri(*i+i»  +Bjh'2[fi(ii+,)l}  . (19) 

The  Wronskian,  W,  is  an  imaginary  constant  (see  eq  (85))  and  is  the  contribution  of  eq  ( 1 7)  at 
the  surface: 

W * -1.45749544104i. 


7.  Some  Effect*  of  Velocity  Structure  on  Low-Frequency  Propagation  in  Shallow  Water,  by  AO  Williams; 
J Acoutt  Soc  Am,  vol  32,  p 363-365,  March  I960. 
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The  depth  functions  are  normalized  by  the  relationship 


Un(zo)Un(z)  = Fn(zO>Fn(z)/Dn-  <20> 

The  functions  F and  F'  used  in  computing  Dn  are  conveniently  assembled  from  the 
elements  of  the  determinant  and  the  coefficients  A,  and  Bj.  This  requires  care  in  developing 
the  computer  code,  because  F is  always  multiplied  by  p and  F'  has  the  term  aj  in  it.  The  sur- 
face differs  from  the  other  layers  in  that  F j is  zero  there  and  F j , by  eq  ( 1 9),  is  a j W.  How- 
ever, because  p\  appears  as  a factor  in  the  coefficients  of  F j , the  actual  value  of  F j at  the 
surface  in  the  computation  is  pj  ajW.  This  factor  of  pj  together  with  the  p^2  needed  for 
orthogonality,  when  squared,  gives  the  p-j  of  eq  (18). 

DETERMINANT 

Normal  modes  are  determined  by  finding  the  eigenvalues  of  a characteristic  equation 
which,  in  turn,  is  obtained  by  setting  a determinant  to  zero.  The  determinant  is  obtained 
from  the  coefficient  matrix  of  a set  of  linear,  homogeneous  equations  expressing  the  bound- 
ary conditions  as  given  by  eq  (10)  - (13).  Since  the  method  of  handling  this  determinant  is 
a central  feature  of  this  normal-mode  program,  it  is  given  in  detail  here. 

The  first  line  of  the  matrix  is  taken  from  eq  (1 1 ) as 

B,  P,  h2  lf,(0)l  + A,  p,  hj  [fj  (0)]  = 0.  (21) 

At  each  profile  interface,  i,  where  i numbers  the  interfaces  below  the  surface  from  1 to 
N-l,  the  two  boundary  conditions  given  by  eq  (12)  and  (13)  are 

BiPjh2  [fj(zi+j)J  +AjPjhj  [Jj(zi+j)]  -Bj+]  pi+1  h2  [?i+j  (Zj+j)l 

-Ai+|  p\+\  hl  Ifi+1  (zi+l)J  = 0 (22) 

and 

Bi  ai  h2  tfi  <zi+l>J  + Ai  ai  hl  l?i(zi+l>]  " Bi+,  ai+1  h'2  [fi+1  (zi+1)) 

“Ai+1  ai+l  hl  tfi+1  ^zi+l^  = 0 • ^23) 

The  coefficients  of  Aj  in  the  first  equation  and  Bj+i  in  the  second  will  be  the  diago- 
nal elements  of  the  matrix.  The  nonzero  elements  of  the  matrix  will  therefore  be  no  more 
than  two  places  from  the  diagonal.  The  matrix  can  be  stored  in  the  computer  in  an  array  of 
size  (2M-1 ) X 4,  where  M is  the  maximum  number  of  layers  in  the  sound-speed  profile.  In 
the  final  layer,  A^  hj  is  omitted,  as  in  eq  (10).  In  the  program,  the  real  and  imaginary  parts 
are  stored  in  separate  arrays. 

The  sparseness  of  the  matrix  permits  efficient  evaluation  by  a triangularization 
process  of  row  reduction.  For  each  pair  of  rows  representing  a pair  of  equations  given  by 
eq  (22)  and  (23),  the  first  element  from  the  first  equation  and  the  first  two  from  the  second 
equation  must  be  set  to  zero  by  subtracting  the  proper  multiple  of  preceding  rows.  The 
determinant  is  then  the  product  of  the  diagonal  elements  of  the  triangularized  matrix.  The 
value  of  the  determinant,  G,  is  used  in  eq  ( 1 5)  to  find  the  roots  by  iteration. 


Note  that  a value  of  v that  makes  this  determinant  zero,  or  near  zero,  ordinarily  is 
zero  because  only  one  diagonal  element  is  very  small.  For  trapped  modes  this  element  is  at 
the  row  representing  the  first  interface  below  the  mode,  ie  the  interface  just  below  the  layer 
of  positive  gradient  in  which  the  sound  speed  is  equal  to  the  mode  phase  velocity.  For 
unstrapped  modes  it  is  usually  the  final  diagonal  element  that  is  small.  Thus  the  layers  in 
which  the  sound  speed  is  greater  than  the  phase  velocity  of  a mode  do  not  greatly  affect  the 
eigenvalue.  Eigenvalues  are  determined  mainly  by  those  parts  of  the  sound-speed  profile 
that  are  less  than  the  phase  velocity. 

When  an  eigenvalue  is  found,  the  coefficients  A,  and  Bj  must  next  be  evaluated.  As 
mentioned  earlier,  one  coefficient  can  be  arbitrarily  chosen.  This  is  done,  and  eq  (21)  is  satis- 
fied by  letting 

Aj  = P|h2lf  j(0)l 
and 

B,  = -P,h,  (f  ,(0)] . (24) 

The  factor  p | is  used  simply  because  the  number  containing  it  is  easily  available  in  the  pro- 
gram. It  is  divided  out  by  the  normalizing  factor,  D.  Eq  (22)  and  (23)  can  then  be  used  to 
evaluate  the  remaining  coefficients,  but  the  triangularized  form  of  the  matrix  yields  the  coef- 
ficients with  less  computation.  If  gu  is  the  element  in  the  ith  row  and  )th  column  of  the  tri- 
angularized matrix,  then  by  Cramers  rule, 

Bi  = Ai-1  ®2i— 2, 2i— 2 *2i-l,2i/Ei 
and 

Ai  = -Ai_i  g2i— 2, 2i— 2 821-1,21-1^’ 

where 

Ei  = «2i-2,2i-l  ®2i— 1 , 2i ” *2i-2, 2i ®2i— 1 , 2i—  1 * (25) 

A simpler  form  is  used  for  B^  in  the  final  layer  since  there  is  no  Ajj  there. 

In  certain  situations  numerical  problems  can  arise  in  evaluating  the  determinant. 

These  require  some  extra  tests  in  the  subroutine  that  makes  the  evaluation.  The  extra  tests 
will  be  discussed  in  the  section,  NUMERICAL  BREAKDOWN.  A more  routine  problem  is 
the  loss  of  accuracy  that  can  arise  in  subtractions  in  the  row  reduction  of  the  matrix.  This 
loss  results  in  less  sharpness  of  convergence  to  a root.  The  size  of  the  determinant.  G,  can  be 
14  orders  of  magnitude  less  at  a root  than  at  the  general  background  near  the  root.  This 
variation  occurs  because  the  modified  Hankel  functions  can  be  computed  to  about  14-place 
accuracy  in  a computer  with  18  decimal  places  available.  Modes  usually  converge  to  10  or 
1 2 places;  thus  a few  places  are  lost  in  evaluating  the  determinant.  In  some  profiles,  usually 
those  with  multiple  ducts  or  those  in  which  propagation  through  bottom  sediments  plays  a 
large  part,  the  convergence  can  be  much  poorer.  Modes  need  to  converge  to  about  4 places 
to  be  reliable  for  computing  losses,  and  convergence  occasionally  fails  to  meet  this  require- 
ment. The  only  current  cure  for  this  loss  in  accuracy  is  to  go  to  higher-precision  arithmetic 
or  to  compute  the  modified  Hankel  functions  to  greater  accuracy.  For  instance,  a standard 
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matrix  triangularization  routine  that  uses  full  row  and  column  pivoting  has  been  tried  with 
no  resultant  increase  in  accuracy. 


FINDING  EIGENVALUES 

There  are  versions  of  this  program  under  development  that  will  locate  the  eigenvalues 
and  do  the  entire  computation  without  user  intervention.  Currently,  however,  these  versions 
are  reliable  only  for  the  simpler  types  of  profiles  - usually  those  with  only  one  duct  - and 
are  not  ready  to  be  reported.  Locating  eigenvalues  with  the  standard  version  of  the  program 
is  discussed  here. 

The  standard  version  of  the  program  requires  the  user  to  find  the  eigenvalues.  In  this 
version,  each  time  an  eigenvalue  is  determined  by  iteration,  the  resulting  value  is  stored  and 
counted  as  an  eigenvalue.  Therefore,  the  user  must  ensure  that  all  iterations  result  in  good 
roots,  that  all  required  modes  have  been  determined,  and  that  no  modes  are  present  more 
than  once.  In  most  cases  the  user  must  expect  to  make  more  than  one  computer  run  to 
obtain  this  result. 

CONTROL  CARDS 

The  user  controls  the  eigenvalue  determination  by  using  any  of  four  different  types 
of  control  cards.  The  first  type  specifies  an  initial  value  for  v and  an  initial  step  size,  Av. 
These  are  both  complex  numbers  with  a real  and  an  imaginary  part.  G is  then  evaluated  at  v 
and  at  v + Av  to  start  the  iteration.  These  are  essentially  the  vj  and  vj+|  of  eq  (15).  If  these 
two  trial  eigenvalues  are  in  the  vicinity  of  a root,  the  iteration  will  converge  to  that  root. 

The  second  type  of  card  specifies  a line  segment  in  the  complex  plane,  along  which  a 
search  for  eigenvalues,  v,  is  made.  The  end  points  of  the  line  arc  given  along  with  the  number 
of  equally  spaced  points  at  which  the  line  is  to  be  divided.  G is  then  evaluated  at  each  suc- 
cessive division  point  along  the  line  until  a relative  minimum  in  IG^I  is  found,  indicating  that 
a root  is  nearby.  The  iterative  process  is  applied  to  find  the  root.  The  initial  step  size,  Av,  is 
first  computed  to  bring  the  second  evaluation  at  v + Av  as  close  as  possible  to  the  true  root. 
This  is  done  by  using  the  point  which  resulted  in  minimum  |G2|  and  the  points  on  either  side 
of  it  to  determine  the  minimum  of  the  parabola  passing  through  them.  If  v - h,  v,  and  v + h 
are  the  three  points  at  which  G was  evaluated,  it  follows  that  the  distance  from  v to  the  mini- 
mum of  the  parabola 

Av  = h(G(v  + h)  - G(v  - h)|  /2l  2G(v)  - G(v  + h)  - G(v  - h)] . (26) 

When  the  iteration  is  complete,  the  eigenvalue  is  recorded  and  the  program  continues 
to  step  along  in  the  direction  of  the  given  line,  checking  again  for  a minimum.  However,  the 
stepping  is  resumed  from  the  newly  located  root  rather  than  from  the  approximate  location 
where  the  minimum  was  detected.  With  this  correction  in  position,  the  designated  line  does 
not  have  to  hug  the  curve  on  which  the  eigenvalues  arc  located  because  it  is  corrected  at  each 
eigenvalue. 

This  method  of  finding  eigenvalues  has  proven  very  successful.  Its  main  utility  arises, 
though,  because  the  eigenvalues  of  the  trapped  modes  have  negligible  imaginary  parts  and  the 
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search  can  be  made  along  the  real  line.  In  simple  profiles  this  can  often  give  a successful  set 
of  modes  on  the  first  try.  Usually,  only  the  three  initial  eigenvalues  need  to  be  located  by 
this  means  because  further  eigenvalues  can  be  located  by  extrapolation  on  the  previous  three. 
This  is  the  function  of  the  third  type  of  control  card. 

The  third  type  of  card  specifies  the  number  of  additional  modes  to  be  determined  by 
extrapolation.  The  starting  value  of  each  eigenvalue  is  determined  by  extrapolating  from  the 
three  most  recently  determined  eigenvalues  to  find  v.  The  step  size,  Av,  is  chosen  as  0.0001 
times  the  distance  between  the  last  two  eigenvalues.  The  exact  eigenvalue  is  then  determined 
by  iteration.  The  extrapolation  is  the  simple  parabolic  form  for  equal  steps: 

v = 3vn-3vn-l  +vn-2-  (27> 

This  method  of  locating  modes  works  well  when  the  modes  lie  along  a smooth  curve, 
as  usually  occurs  for  single  ducts.  But  this  relationship  does  not  always  occur  for  profiles 
with  multiple  ducts. 

The  final  control  card  is  punched  by  the  program  when  requested  and  contains  the 
correct  eigenvalue  to  full  precision.  Upon  encountering  this  card,  the  program  does  not 
iterate,  but  instead  evaluates  G for  this  eigenvalue  and  stores  this  value  of  G as  the  next 
eigenvalue.  A deck  of  such  completed  eigenvalues  can  be  stored,  saving  the  expense  of 
recomputing  the  eigenvalues  for  a given  profile  and  frequency. 

ITERATION  TERMINATION 

A full  description  of  the  iteration  of  eq  ( 1 5)  should  include  the  method  of  termina- 
tion. The  usual  criterion  for  stopping  is  that  G fails  to  become  smaller.  As  G approaches 
minimum  size,  however,  round-off  error  can  act  as  noise  so  that  G is  no  longer  a predictable 
function  of  v.  The  denominator  of  eq  (15)  can  then  be  very  small  by  chance,  resulting  in  a 
large  value  for  6j.  If  this  happens,  the  next  value  of  v,  which  was  as  near  to  the  root  as  pos- 
sible, will  be  far  away.  A much  better  convergence  criterion  is  that  6j  has  reached  a minimum 
in  absolute  value.  In  the  program,  iteration  is  stopped  when  |62|  exceeds  the  previous  value 
by  a factor  of  2.  However,  this  criterion  is  not  applied  until  three  iterative  steps  have  been 
completed,  to  permit  the  process  to  become  well  established.  An  upper  limit  of  15  iterative 
steps  is  permitted.  We  have  not  found  an  improvement  on  the  root  after  1 5 steps. 


* 

I 
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SOUND  SPEED  PROFILE 

The  normal  mode  program  requires  as  inputs  the  depth  of  each  layer  and  the  sound 
speed  and  sound  speed  gradient  at  the  top  of  each  layer.  These  variables  are  mapped  into  the 
dimensionless  internal  variables  of  the  program  by  eq  (7).  The  purpose  of  the  sound  speed 
profile  processing  portion  of  the  program  is  to  accept  the  profile  parameters  in  a form  con- 
venient for  the  user  and  to  translate  them  into  the  required  sound  speeds  and  gradients. 

The  first  function  of  the  processing  program  is  to  make  the  sound  speed  continuous 
at  interfaces.  This  is  done  simply  by  using  the  sound  speed  at  the  bottom  of  one  layer  as  the 
sound  speed  at  the  top  of  the  next.  It  may  be  necessary  to  compute  the  sound  speed  at  the 
bottom  of  the  layer.  The  necessary  parameters  will  have  been  given.  Occasionally  a 
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discontinuity  in  sound  speed  is  required,  as  when  modeling  an  interface  between  water  and 
sediment.  The  user  indicates  this  by  specifying  the  sound  speed  at  the  top  of  the  layer.  If 
left  blank,  the  program  provides  the  sound  speed  necessary  for  continuity. 


A second  function  of  the  processing  program  is  to  permit  a layer  to  be  defined  by  the 
sound  speed  at  top  and  bottom  of  the  layer  rather  than  by  one  sound  speed  and  one  gradient. 
Note  that  the  profile  form  as  given  by  eq  (5)  is  a two-parameter  curve. 

The  last  layer  extends  to  infinite  depth,  so  a gradient  must  be  specified  at  the  top  of 
it.  However,  this  gradient  can  be  specified  by  giving  a depth  and  sound  speed  point  below  the 
last  layer.  The  program  handles  this  by  checking  to  see  if  the  gradient  of  the  last  given  layer 
is  unspecified.  If  it  is,  the  number  of  layers  is  reduced  by  one,  which  causes  the  last  layer  to 
be  only  the  required  extra  point  determining  the  final  gradient.  This  final  gradient  must 
always  be  negative,  as  is  required  by  the  boundary  conditions.  The  program  user  must  ensure 
that  this  gradient  is  negative  and  that  no  gradient  is  zero.  A zero  gradient  will  appear  in  the 
denominator  of  eq  (7). 

These  functions  of  the  profile  processing  program  are  relatively  simple,  but  an  addi- 
tional capability  used  to  model  sediment  bottoms  greatly  increases  the  complexity  of  the 
program.  The  capability  required  is  to  specify  the  absorption  in  a layer  by  adding  an  imagi- 
nary part  to  the  sound  speed.  In  older  versions  of  this  normal  mode  program  an  imaginary 
part,  expressed  as  an  absorption  coefficient,  could  be  added  to  the  sound  speed  at  the  top  of 
the  layer.  This  imaginary  part  is  small  compared  to  the  real  part.  Since  the  gradient  was 
assumed  real  at  the  top  of  the  layer,  the  imaginary  part  was  initially  not  changing  with  depth 
and  it  usually  changed  only  a minor  amount  through  the  depth  of  the  layer.  However,  this 
small  change  could  not  always  be  relied  upon.  Aiso  Hamilton  (ref  8)  has  published  data  on 
absorption  gradients  in  sediment  layers,  so  more  precise  control  of  this  part  of  the  sound 
speed  function  is  needed  to  model  sediment  layers.  Therefore,  a more  comprehensive  profile 
processing  routine  has  been  incorporated  in  the  normal  mode  program.  This  curve-fitting 
process  is  described  below. 

The  following  quantities  can  be  input  for  each  layer  depth  starting  at  the  surface: 

Depth  of  top  of  the  layer 

Sound  speed  at  top  of  layer 

Sound  speed  at  bottom  of  layer 

Real  part  of  sound  speed  gradient  at  top  of  layer 

Attenuation  in  loss  per  km  at  the  top  of  the  layer 

A similar  attenuation  at  the  bottom  of  the  layer 

Density  in  the  layer 

The  density  is  a constant  in  the  layer  and  as  such  requires  no  further  curve  fitting.  Redun- 
dant parameters  are  left  blank  on  input  cards.  In  some  cases  negative  values  serve  as  flags  to 
indicate  specific  treatment.  For  instance  a negative  value  of  absorption  at  the  top  of  a layer 


8.  Sound  Attenuation  as  a Function  of  Depth  in  the  Sea  Floor,  by  EL  Hamilton;  J Acoust  Soc  Am,  vol  59, 
p 528-535,  March  1976. 
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directs  the  program  to  use  the  same  imaginary  part  of  sound  speed  as  occurred  at  the  bottom 
of  the  previous  layer.  Similar  flags  at  the  bottom  of  a layer  are  discussed  later. 


Absorption  per  Hz  is  given  in  units  of  decibels  per  km  (or  kiloyard).  The  quotient  of 
absorption  over  frequency  is  used  because  Hamilton  (ref  8)  usually  considers  absorption  (or 
attenuation)  as  proportional  to  frequency  with  a coefficient  k.  We  use  the  symbol  h instead. 
That  is. 


or  = hf. 

We  interpret  a to  be  in  units  of  dB  per  km  and  f in  Hz,  whereas  Hamilton  uses  dB  per  m and 
kHz;  but  the  coefficients  h and  k remain  equal. 

The  complex  wave  number  in  layer  i is  represented  as 

kj  = cj/q 

= uReC|Cr2-iwlmC|C|"2.  (28) 

A plane  wave  will  be  attenuated  a dB  per  km  if 
Imkj  = -a/(20000  log  e) 

= -irAf,  (29) 

where 


A = h/(20  000  tr  log  e). 

By  equating  the  imaginary  part  of  kj  in  eq  (28)  and  (29),  the  imaginary  part  of  Cj  is  found  to 
be  as  follows: 

ImCi=l/A-[l/A2-(ReCi)2lW  (30) 

If  a is  zero,  which  is  the  case  usually  used  in  water  layers,  eq  (30)  cannot  be  used;  but  the 
imaginary  part  of  C is  then  simply  zero.  These  two  cases  are  treated  separately  in  the  program. 

When  sound  speed  is  given  at  the  top  and  bottom  of  layer  i,  the  imaginary  parts  of 
the  sound  speeds  are  determined  by  eq  (30)  and  the  only  curve  fitting  task  is  to  determine  the 
gradient  yj.  Solving  eq  (5)  for  yj, 

Y,  = C,(C*,  - Cj2)/2C2  jtx,,., (31) 

The  gradient  is  a complex  number  since  the  C’s  here  are  complex.  The  z’s  are  real. 

A second  version  of  this  computation  arises  if  the  gradient  is  required  to  be  a real 
number.  In  this  case,  which  is  used  to  match  older  versions  of  the  program,  an  additional 
parameter  must  be  left  unspecified  and  this  parameter  is  Im  Cj+j.  This  is  equivalent  to  having 
the  sound  absorption  at  the  bottom  of  the  layer  unspecified.  Therefore,  a negative  number 
input  for  this  parameter  is  used  as  a flag  to  call  for  this  particular  fitting  procedure. 
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For  this  situation,  given  Re  Cj,  Im  Cj,  Re  Cj+  j , and  making  7j  real,  the  determination 
of  and  Im  Cj+j  is  not  simple.  When  7,  is  eliminated  from  the  real  and  imaginary  parts  of 
eq  (31),  a quartic  equation  in  Im  Cj+j  results.  Rather  than  derive  an  algebraic  solution  to 
this  equation,  it  is  solved  by  iteration  under  Newton’s  method.  A good  first  guess  at  the 
solution  is  Im  Cj+i  2 Im  Cj.  Four  iterations  usually  give  an  accurate  root.  The  equation  is 

Im  Cj  (Im  Cj+j)4  + [lm(Cj)^  + 2(Re  Cj+j)2  Im  CjJ  (Im  Cj+j)2 
+ 2Re  Cj+j  Re(Cj)^  Im  Cj+j 

+ Im  Cj(ReCi+1  )4  - (Re  Ci+, )2  Im(C3)  = f(Im  Ci+1 ).  (32) 

The  root  is  then  found : 

(ImCj+1).  = (ImCi+1)._1  -f If. 

The  gradient,  7,  is  next  given  by  the  relationship 

7j  = Im  Cj  [(Re  Ci+1  )2  - (Im  Ci+1)2]  + 2Re  Cj  Re  Ci+1  Im  Ci+1  - Im(C3) 

[4  Re  Ci+1  Im  Ci+ j (Zj+,  - Zj)] . (33) 

Because  the  root  of  eq  (32)  may  not  be  exact,  Im  75  may  not  be  exactly  zero.  This  slight 
error  can  be  transferred  to  Cj+j  by  using  the  computed  real  7j  to  recompute  Cj+j.  This  is 
done  in  the  program  by  transferring  to  a portion  of  the  program  already  designed  to  do  this. 

When  sound  speed  and  gradient  at  the  top  of  the  layer  are  given,  the  parameters 
required  by  the  program  are  all  given.  The  sound  speed  at  the  bottom  of  the  layer  is  rou- 
tinely computed,  however,  because  it  may  be  required  to  make  the  next  layer  continuous. 
Equation  (5)  is  used  to  determine  the  sound  speed  at  depth  zj+j,  which  is  the  depth  of  the 
bottom  of  the  layer.  This  is  straightforward,  but  several  complications  arise.  Only  the  real 
part  of  the  gradient  at  the  top  of  the  layer  is  used  as  an  input  because  situations  have  not 
arisen  that  require  that  the  imaginary  part  of  the  gradient  be  specified.  Often  the  attenua- 
tion is  given  at  both  top  and  bottom  of  the  layer.  That  is,  Re  Cj,  Im  Cj  and  Re  7j  are  given, 
plus  a relationship  between  Re  Cj+j  and  Im  Cj+j.  The  imaginary  part  of  the  gradient,  Im  74, 
must  be  determined  as  well  as  both  real  and  imaginary  parts  of  the  sound  speed  at  the  layer 
bottom.  The  derivation  of  this  case  is  not  trivial. 

One  relationship  between  the  real  and  imaginary  parts  of  the  sound  speed  is  given  by 
eq  (28)  and  (29).  From  these  equations  at  Cj+j  we  derive 

A(T-i)  = 2/Cj+j,  (34) 

where 

T=  ReCi+j/ImCj+j. 

Substituting  this  expression  for  Cj+j  into  eq  (31)  and  equating  real  parts  gives  a quadratic 
expression  for  T which  has  a usable  root  of 
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(35) 


Re(C3)T  * -Im(C3)  - 


(Im(C3)] 2 + Re(C?)B 


wnere 


B * Re(C3)  - 8 Re  7j(zi+i  - Zj)/A2  + 4 Re  Cj/A2. 

From  eq  (34), 

ReCi+1  = 2T/A(T2  + I) 
and 


= RCj^.j/T. 


(36) 


The  gradient  can  now  be  evaluated  by  eq  (31)  to  find  its  imaginary  part. 


Equations  (34)  and  (35)  cannot  be  used  if  the  attenuation  at  the  bottom  of  the  layer 
is  given  as  zero.  Therefore  an  alternate  form  must  be  used.  This  form  is  much  simpler  than 
the  previous  case,  since  Cj+|  is  real. 

* 


Ci+1  = 


,3 


Re(Cj ) / (Re  Cj  - 2 Re  (zj+,  - Zj)l 


Im  7j  *=  ( lm  Cj  - lm(Cj  )/Cj+1 ) [2(zi+,  - zj)] 


-1 


(37) 

(38) 


Finally,  if  the  special  case,  7j  real,  is  specified  by  inputting  a negative  value  for 
absorption,  eq  (31)  can  be  used  directly  to  give 

Ci>l  "CjVlCj- 27j(*j+| 


(39) 


To  evaluate  the  square  root,  let 


- a + bi. 


Then 


and 

lm  Cj+i  = b/2  Re  Cj+|. 


(40) 

(41) 


NUMERICAL  BREAKDOWN 

A situation  arises  frequently  in  which  a very  small  depth  function  must  be  computed 
from  the  difference  of  two  large  numbers.  A wrong  answer  results  if  this  accuracy  loss 
exceeds  the  word  size  of  the  computer.  The  best  way  that  has  been  found  to  avoid  this  is  to 
check  for  it  within  the  program  and  arbitrarily  replace  the  wrong  number.  In  checking  for 
this,  a constant,  called  T-lim  in  the  computer  program,  is  compared  to  the  argument  of  th' 
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modified  Hankel  functions  or  to  the  argument  of  the  exponential  function  within  modified 
Hankel  functions.  A T-lim  value  of  25.0  is  used  in  the  program,  but  a smaller  number  occa- 
sionally is  required.  The  program  user  can  alter  T-lim  by  appropriate  input  cards  (Key  8 = 1 
followed  by  a new  value  of  T-lim).  The  next  few  paragraphs  demonstrate  the  symptoms  of 
this  problem,  so  as  to  assist  a user  in  recognizing  the  problem.  The  remainder  of  this  section 
describes  the  modifications  that  have  been  made  to  the  computer  program  to  correct  this 
loss  of  accuracy. 

The  solid  line  of  figure  1 shows  a simple  surface  duct  and  the  phase  velocities  of  the 
first  three  modes  at  3 kHz.  For  this  profile,  the  depth  function  of  mode  1 is  shown  in  fig- 
ure 2.  The  solid  line  is  the  depth  function  as  computed  by  a program  that  does  not  correct 
for  numerical  breakdown.  The  dashed  line  shows  the  correct  depth  function  below  a depth 
of  71  m.  This  result  was  determined  from  Airy  functions,  not  from  the  program.  Between 
depths  of  71  to  100  m,  the  program  cannot  compute  the  depth  function  accurately.  In  the 
second  layer,  which  starts  at  a depth  of  100  m,  the  function  can  be  computed  accurately  but 
it  is  incorrectly  placed  by  the  boundary  condition  that  requires  the  depth  functions  to  be 
continuous  at  interfaces.  The  slope  of  the  depth  function  was  correctly  computed  as  indi- 
cated by  the  identical  shape  of  the  three  depth  functions  in  the  second  layer.  The  shape  is 
such  as  to  make  the  correct  depth  function  continuous  in  slope  across  the  interface. 

The  breakdown  in  accuracy  at  a depth  of  71  m occurred  when  f had  a value  of -8.4. 
(?  is  given  by  eq  (7)  and  is  the  argument  of  the  modified  Hankel  functions.)  A negative  value 
of  ? occurs  when  the  mode  phase  velocity  is  less  than  the  speed  of  sound.  Since  the  ray  of 
the  same  phase  velocity  cannot  reach  such  a region,  the  sound  field  there  is  a diffracted  field. 
The  mode  depth  function  is  therefore  small  at  such  depths.  In  the  figure,  the  depth  func- 
tion amplitude  at  the  breakdown  point  is  about  7 orders  of  magnitude  (or  in  terms  of  propa- 
gation loss,  140  dB)  down  from  its  maximum.  Equations  (62),  (66),  and  (68),  which  will  be 
given  for  the  modified  Hankel  functions,  indicate  that  the  argument  of  the  exponential  term 
is  2/3(8. 4)3/2,  or  16.2.  The  functions  hj  and  h2  will  thus  be  about  107  in  magnitude  at  a 
depth  of  71  m.  These  large  values  and  their  small  difference  account  for  the  approximate 
accuracy  loss  of  14  decimal  places,  which  is  the  general  accuracy  of  the  modified  Hankel 
functions. 

Incorrect  behavior  in  the  depth  function  usually  occurs  when  ? is  about  -8.4.  In 
some  more  complicated  profiles,  however,  where  accuracy  is  lost  in  row  reduction  of  the 
determinant,  the  depth  functions  may  become  incorrect  at  values  of  ? that  are  less  in  abso- 
lute value.  When  this  problem  occurs  it  can  be  diagnosed  by  plotting  the  depth  function  of 
the  mode  and  noting  the  steep  positive  slope  through  some  depth  interval  as  in  figure  2. 

When  that  occurs,  the  value  of  T-lim  should  be  decreased. 

Incorrect  depth  functions  can  cause  errors  in  propagation  loss  computations  in  two 
ways.  In  figure  2,  the  solid-line  depth  function,  because  of  its  large  size,  can  cause  losses  to 
be  too  low  at  a depth  of  around  1 00  m.  The  second  error  would  occur  if  the  duct  were 
deeper,  say  1 10  m.  At  this  depth  the  erroneous  segment  of  depth  function  in  figure  2 would 
reach  a value  of  about  10“I , where  it  would  be  larger  than  the  correct  lobe  of  the  depth 
function  near  the  surface.  With  this  extra  area  under  the  curve,  the  normalizing  factor  would 
be  increased  significantly  and  would  reduce  the  size  of  this  entire  depth  function.  Thus, 
losses  near  the  surface  would  be  larger  because  of  the  loss  in  size  of  mode  1 . 
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Figure  1 . A two-layer  sound  speed  profile  for  a 
surface  duct.  The  phase  velocities  of  the  first 
three  modes  at  3 kHz  are  marked.  The  broken 
line  shows  a modification  of  the  upper  layer  to 
prevent  numerical  breakdown  in  mode  I . 


The  standard  correction  to  mode  1 is  shown  in  figure  2 by  the  dot-dashed  line.  In  the 
depth  interval  where  f < -8.4,  the  function  is  set  to  zero.  The  values  of  the  depth  function 
at  greater  depths  result  from  a modification  in  the  values  used  in  the  determinant. 

The  corrected  values  in  figure  2 are  not  equal  to  the  true  value  of  the  depth  function, 
but  they  are  small  enough  that  they  do  not  alter  the  propagation  loss  to  a tenth  of  a decibel 
when  a full  set  of  modes  is  used.  When  the  source  or  receiver  is  at  a depth  where  such  cor- 
rections are  necessary,  the  mode  can  be  omitted  from  the  computation.  Thus,  properly 
omitting  modes  would  solve  the  above  problems  except  for  the  cases  where  the  normalizing 
factor,  D,  is  affected.  In  these  cases,  losses  cannot  be  computed  accurately  without  the 
corrections. 

PROGRAM  MODIFICATIONS 

The  modification  is  approximately  equivalent  to  modifying  the  sound  speed  profile 
as  shown  by  the  broken  line  in  figure  1 . In  effect,  the  sound  speed  is  not  allowed  to  become 
enough  greater  than  the  phase  velocity  of  the  mode  being  considered  to  cause  problems. 

Tile  limi’ation  on  £ is  accomplished  at  three  different  places  in  the  normal  mode  pro- 
gram. It  is  not  clear  that  this  is  the  best  way  to  handle  the  problem  and  it  may  be  redundant, 
but  it  appears  to  be  an  adequate  solution.  These  three  corrections  will  be  described  next. 
Finally  a correction  to  the  determinant  program  is  described  which  is  necessary  because  the 
limiting  of  f can  cause  false  zeroes  in  the  determinant. 

In  the  subroutine  SETUP  the  elements  of  the  determinant  are  computed  by  deter- 
mining f at  the  top  and  bottom  of  each  layer  and  then  calling  the  modified  Hankel  function 
program.  At  the  top  ot  each  layer.  Re  f is  set  to  -7.5  if  its  value  was  less.  However,  this  is 
done  in  an  iterative  loop  in  which  the  real  part  of  w/Cj  in  eq  (7)  takes  on  a sufficiently  larger 
value  while  its  imaginary  part  is  fixed.  This  is  done  to  retain  the  absorptive  properties  of  a 
layer  when  its  sound  speed  is  in  effect  being  reduced.  It  has  been  found  unnecessary  to  make 
the  above  constant,  -7.5,  a t unction  ot  T-lini  which  the  user  can  vary,  because  an  oversized 
value  at  the  top  of  a layer  is  not  as  critical  as  at  the  bottom.  At  the  bottom  of  a layer,  several 
tests  are  made.  If  the  real  part  of  £ has  decreased  past  the  limit  at  some  depth  between  the 
top  and  bottom  of  the  layer,  it  is  set  at  the  limit.  This  limit,  called  S-lim  in  the  program,  is 
related  to  T-lim  by  the  relationship 

S“-(T)2/3  (42) 

where  S and  T are  the  two  limits.  If  Re  f is  less  than  -7.5  throughout  the  layer,  it  is  simply 
set  at  -7.5.  Such  a layer  has  negligible  effect  on  a mode. 

In  program  MAIN  at  the  location  where  depth  functions  are  computed  for  given 
depths,  a process  similar  to  that  above  is  used.  To  evaluate  the  depth  function  in  a given 
layer,  £ is  first  evaluated  at  the  top  of  the  layer.  The  real  part  of  f is  theq  limited  as  in  the 
program  SETUP  above.  Next  f is  evaluated  at  the  given  depth  by  adding  the  depth-dependent 
part  onto  the  value  at  the  top  of  the  layer  which  may  be  the  modified  value.  If  this  final 
value  is  less  than  S-lim,  the  depth  function  is  set  to  zero.  If  it  is  greater  than  S-lim.  the  func- 
tion is  computed  in  the  usual  way. 
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The  imaginary  part  of  f can  be  large  if  the  eigenvalue  has  a large  imaginary  part  or  if 
the  speed  of  sound  in  the  layer  has  a large  imaginary  part.  When  this  happens  the  imaginary 
part  of  2/3  f3/2t  which  appears  as  an  exponential  in  the  modified  Hankel  functions,  may 
become  large  in  absolute  value  even  though  Re  f has  been  limited.  A final  check  is  therefore 
made  before  the  exponential  is  computed.  If  Im  f^/2  js  greater  than  T-lim,  f is  reduced  in 
amplitude  to  the  size  at  which  it  will  equal  T-lim.  The  angle  of  J in  the  complex  plane  is 
preserved. 

This  limitation  of  the  exponential  can  be  viewed  in  another  way.  In  a following  sec- 
tion the  two  components  of  the  modified  Hankel  functions,  Fj  and  F2,  eq  (68)  and  (69), 
have  exponential  terms  whose  arguments  are  equal  and  opposite  in  sign.  When  these  argu- 
ments have  magnitude  of  2/3  T-lim,  they  differ  in  size  by  1 5 decimal  places,  which  is  near 
the  1 8-decimal-place  word  size  of  the  machine.  The  ability  to  compute  the  difference  in 
these  two  terms  is  essentially  the  same  as  the  ability  to  compute  the  depth  function  accurately. 

PREVENTING  ZEROES  IN  THE  DETERMINANT 

Placing  limits  on  f can  cause  problems  in  the  determinant  because  f may  be  set  equal 
to  S-lim  at  several  interfaces.  The  equations  that  arise  for  matching  boundary  conditions 
may  then  be  identical  for  these  interfaces  and  may  therefore  fail  to  be  linearly  independent. 
The  triangularized  determinant  will  thus  have  zeroes  on  the  diagonal  at  positions  equivalent 
to  interfaces  that  do  not  have  real  physical  importance  for  the  mode.  These  will  prevent 
location  of  the  significant  “zeroes”  or  roots.  These  artificial  zeroes  must  be  removed. 

The  artificial  zeroes  are  detected  and  removed  in  the  subroutine  DET,  which  evaluates 
the  determinant.  If  four  elements  from  the  matrix  have  the  configuration 

a b 


c d 


and  c is  to  be  set  to  zero  by  row  reduction,  d will  be  replaced  by  a value,  x,  as  follows: 
x = d - bc/a. 

If  d is  located  on  the  diagonal,  complete  loss  of  accuracy  is  checked  for  by  computing 
s=|x2|/|d2|. 

If  s is  less  than  10“^,  x is  not  used;  instead,  d is  replaced  by  10~*  ^d.  Note  that  this  substi- 
tution will  occur  when  x is  zero,  thus  preventing  zeroes  on  the  diagonal.  The  power  of  ten, 
-1 7,  is  chosen  to  be  near  the  total  word  size  of  1 8 decimal  places. 

The  above  substitution  prevents  sudden  jumps  in  the  value  of  the  determinant  when 
all  precision  is  lost  at  one  step  in  the  evaluation.  This  is  important  for  the  mode  search  rou- 
tine which  detects  roots  by  looking  for  minima  in  a series  of  values  of  the  determinant  while 
one  parameter  is  incremented  slowly.  A sudden  jump  will  often  produce  a relative  minimum 
which  will  be  falsely  interpreted  as  a root.  At  true  roots,  one  or  more  elements  along  the 
diagonal  are  small,  but  not  as  small  as  those  checked  for  here. 
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REFLECTION  COEFFICIENTS  AND  OTHER  AUXILIARY  OUTPUTS 


Once  the  depth  functions  of  a mode  have  been  determined,  it  is  relatively  easy  to 
compute  reflection  coefficients  at  any  interface.  Therefore,  a subroutine  called  RCOEF  has 
been  added  to  the  program  which  will  compute  and  print  out  reflection  coefficients  if 
requested  by  the  use  of  control  key  3.  If  key  3 is  set  to  1 , the  reflection  coefficients  at  all 
interfaces  are  computed.  If  set  to  a number,  n,  greater  than  1 , the  coefficient  is  computed 
at  the  nth  interface  only,  where  the  surface  is  the  first  interface. 

The  printout  includes  the  phase  as  well  as  the  amplitude  of  the  reflection  coefficient 
and  the  grazing  angle.  The  grazing  angle,  0,  of  the  equivalent  rays  is  computed  from  the 
mode  phase  velocity  and  the  sound  speed,  c,  at  the  bottom  of  the  layer,  by  Snell’s  law: 

6 - cos-'  (c/v). 

The  grazing  angle  is  computed  only  if  the  phase  velocity  is  greater  than  the  sound  speed  at 
the  interface,  since  otherwise  the  equivalent  ray  does  not  reach  the  interface. 

The  reflection  coefficient  is  derived,  following  Bucker  (ref  9),  by  assuming  that  an 
isospeed  layer  exists  for  a small  depth  just  above  the  interface.  In  this  layer  the  depth  func- 
tion can  be  written  as 

f(z)  = AeUz  + Be-ilz,  (43) 

where  1,  the  vertical  component  of  the  mode  wave  number,  is  given  for  mode  n by 


and 

kj  = «/cbi> 

where  c^j  is  the  sound  speed  at  the  bottom  of  layer  i.  The  derivation  now  consists  of  identi- 
fying A and  B as  the  pressures  of  the  upgoing  and  downgoing  waves  at  the  bottom  of  the 
layer;  thus  the  reflection  coefficient 

R = A/B. 

A and  B are  evaluated  by  making  f and  its  derivative  at  the  interface  between  this  small  iso- 
speed layer  and  the  regular  profile  continuous  with  the  normal  mode  depth  functions.  The 
thickness  of  the  isospeed  layer  is  then  allowed  to  approach  zero,  giving  the  desired  value  of 
R.  If  F and  F'  are  the  normal  mode  function  and  its  depth  derivative  at  the  interface  depth 
defined  by  eq  (9)  and  (19),  the  reflection  coefficient  resulting  from  the  above  derivation  is 
as  follows: 

R = (ilF  + F')/(ilF  - F').  (45) 

This  coefficient  is  a complex  number.  Loss  per  reflection  is  given  by  20  times  the 
log  of  the  absolute  value.  The  phase  gives  the  phase  shift  that  an  equivalent  ray  would 


9.  Sound  Propagation  in  a Gunnel  with  Lossy  Boundaries,  by  HP  Bucker;  J Acoust  Soc  Am,  vol  48, 
p 1 187-1 194,  November  1970. 


20 


experience  upon  reflection.  Figure  3 is  an  example  of  the  use  of  this  computation.  It  shows 
phase  and  amplitude  of  the  reflection  coefficient  in  shallow  water  over  a sandy-silt  sediment 
lying  over  rock.  The  frequency  is  1 500  Hz.  Reflections  are  given  only  at  discrete  points 
determined  by  the  individual  modes. 

The  model  in  figure  3 is  for  a liquid  bottom.  That  is,  no  rigidity  is  supplied  in  this 
program  and  the  sound  speed,  density,  and  attenuation  determine  the  reflection  coefficients. 

The  reflection  coefficients  computed  by  eq  (45)  can  be  closely  approximated  by 
dividing  the  mode  attenuation  by  the  loop  length  of  the  corresponding  ray.  The  loop  length 
must  be  determined  from  ray  theory  for  the  ray  of  the  same  phase  velocity  or  vertexing  veloc- 
ity. However,  an  interesting  analog  of  the  ray  loop  length  is  the  intermode  interference  length. 
This  is  discussed  by  Guthrie  (ref  10).  Specifically, if  the  difference  between  eigenvalues,  ReXj, 
for  two  adjacent  modes  is  AX,  the  interference  length  1 = 2w/AX.  This  distance  will  usually 
equal  the  ray  loop  length  for  some  ray  with  phase  velocity  between  that  of  the  two  modes. 

As  each  mode  after  the  first  is  computed,  the  length,  1,  is  computed  and  printed  out. 
Also  routinely  printed  out  for  each  mode  is  the  mode  damping  or  mode  attenuation  coeffi- 
cient, in  units  of  dB  per  km.  This  attenuation,  aj,  is  computed  from  the  relationship 

oij  = -1000  Im  Xj  log  | q e 

= - 8686  1m  Xj. 

This  quantity  multiplied  by  range  gives  the  damping  of  mode  i,  in  dB. 


10.  The  Connection  Between  Normal  Modes  and  Rays  in  Underwater  Acoustics,  by  K.M  Guthrie;  J of  Sound 
and  Vibration,  vol  32,  no  2,  p 289-293, 1974. 
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Figure  3.  A shallow-water  profile  with  resulting  phase  and  amplitude  of  the  reflection  coefficient  at  1 .5  kHz. 
Parameters  at  the  top  of  the  sediment  layers  are  as  follows:  1st  layer  - c = 1606.45  m/s,  y * l.5s'*,a  = 
0.18dB/m,  p ■ 1.68;  2nd  layer  -c  * 1684.0  m/s,  y = 1.5s"*, a * 1.10dB/ni,p  = 1.91;  final  layer  - y = -0.1. 
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COMPUTATION  OF  THE  MODIFIED  HANKEL  FUNCTIONS 


Most  of  the  computer  time  required  to  determine  eigenvalues  and  compute  depth 
functions  is  spent  in  evaluating  the  modified  Hankel  functions  of  order  1/3.  For  this  reason, 
minimizing  computer  time  in  evaluating  these  functions  is  desirable.  Gaining  as  many  places 
of  accuracy  as  possible  is  even  more  important.  The  average  normal  mode  computation  will 
have  many  modes  that  can  be  determined  to  far  greater  accuracy  than  is  required  to  obtain 
0.1  dB  accuracy  in  the  propagation  loss.  However,  there  are  usually  some  and  often  many 
modes  in  which  many  places  of  accuracy  are  lost  in  evaluating  the  determinant.  Therefore, 
maximum  accuracy  in  the  modified  Hankel  functions  is  required  to  extend  the  range  of 
cases  for  which  computations  can  be  carried  out  successfully. 

Optimization  of  the  program  is  a function  of  the  computer  word  length.  The  pro- 
gram given  in  this  report  is  for  the  UNI  VAC  1110  with  60  bits  word  length  in  double  preci- 
sion or  18.1  decimal  places.  This  section  gives  the  equations  and  computational  techniques 
that  are  required  to  optimize  this  program  for  different  computer  word  lengths.  Complete 
details  of  the  functions  are  given  in  reference  1 1. 

The  Airy  functions  Ai(Z)  and  Bi(Z)  can  be  used  instead  of  the  modified  Hankel 
functions  h j and  h2-  However,  since  h2  is  ideally  suited  to  matching  the  boundary  condi- 
tions at  great  depth  as  formulated  in  this  normal  mode  program,  hj  and  I12  are  used  here. 
The  relationship  between  them  is  as  follows: 


h j (z)  = k [Ai  (— z)  - i Bi  (-z)J 

(46) 

h2(z)  = k*  [Ai  (— z)  + i Bi  (-z)J  , 

(47) 

where 


k = (3/2)^  (1  - iv^/3),  and  k*  is  the  complex  conjugate  of  k. 

In  this  section  z will  be  the  argument  of  the  functions  h j and  h2-  For  small  values 
of  |z|,  h j and  h2  are  computed  by  power  series  expansions.  For  large  values,  an  asymptotic 
expansion  is  used.  In  the  past  the  asymptotic  series  was  expanded  directly.  However,  a 
continued  fraction  expansion  has  been  found  to  give  both  shorter  running  time  and  better 
accuracy. 

Figure  4 shows  a line  in  the  complex  plane  which  divides  the  plane  into  two  parts. 
For  values  of  z within  the  line,  the  power  series  method  is  used.  When  z is  outside  the  line, 
the  continued  fraction  method  is  used.  This  line  is  a function  of  computer  word  length,  and 
the  method  of  determining  it  will  be  given  after  the  two  methods  have  been  treated.  The 
accuracy  of  the  methods  is  also  treated. 

The  program  has  a parameter  called  IH  in  the  FORTRAN  call  statement  which  con- 
trols which  functions  are  computed.  If  IH  is  set  to  zero,  both  functions  and  their  derivatives 
are  computed.  If  IH  is  set  to  1 , only  the  functions  are  computed.  If  it  is  set  to  2,  only  h2 
and  its  derivative  are  computed. 


1 1.  Tables  of  the  Modified  Hankel  Functions  of  Order  One-Third  and  their  Derivatives,  Harvard  University 
Computation  Laboratory;  Harvard  University  Press,  Cambridge  MA,  1945. 
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Figure  4.  Line  in  complex  plane  dividing  the  arguments  for  which  the 
modified  Hankel  functions  are  computed  by  ( 1 ) power  series  (inside) 
and  (2)  asymptotic  expansion  evaluated  by  continued  fractions 
(outside). 


POWER  SERIES  EXPANSION 

In  this  expansion  hj  and  h2  are  given  by  two  auxiliary  functions  f and  g as 


h,  (z)  = g + i(3)-1/2(g-20  (48) 

h2(z)  = g-i(3r1/2(g-2f)  • (49) 

The  auxiliary  functions  are  given  by  the  expressions 
M 

f = A £ a,,,  Xm  (50) 

m=0 

M 

g = BZ  2 bmXm  - (51) 

m=0 


where  X = z^,  A = 2^/^/[r(2/3)]  and  B = 2^/[32^  T (4/3)).  The  derivatives  h'j  (z)  and 
h'2  (z)  can  be  derived  by  straightforward  differentiation  of  eq  (50)  and  (51)  to  give 
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(52) 


f 


g 


M 

-Az2  y cmxm 
m=0 
M 

B J dmX">  . 
m=0 


(53) 


The  coefficients  of  eq  (50)  - (53)  are  given  by  recursion  relations  where  aQ  = 1 , aj  = 1/3!, 
a2  = + l • 4/6!,  a3  = -l  • 4 • 7/9! 

am  = _am-l/(3m)(3m-l)  (54) 


b0  = I,  bj  = -2/4!,  b2  = + 2 • 5/7*,  b3  =-2*5*  8/10! 


bm  = -bm_1/(3m)(3m+l)  (55) 

cq  = 3/3!,  Cj  = -6  • 1 • 4/6! 

cm  = ~cm_i/3m  (3m+2)  (56) 

d0=l,d,  =-4-2/4! 

dm  = -dm_i/3m  (3m-2)  (57) 

It  is  important  for  efficient  computation  that  the  number  of  terms  M be  no  larger 
than  necessary.  In  the  current  program  the  same  value  of  M is  used  in  all  four  sums.  This  is 
done  because  the  optimum  number  never  differs  by  more  than  one  in  the  four  cases  and  the 
determination  by  table  look-up  of  four  M’s  often  would  take  longer  than  computing  any 
unnecessary  terms.  M for  each  series  is  determined  so  that  adding  additional  terms  will  not 
change  the  answer.  Then  the  most  stringent  of  the  four  conditions  is  tabulated  and  used. 

A precise  determination  of  the  number  of  terms  to  use  requires  a knowledge  of  the 
size  of  the  largest  single  term  in  the  sum.  When  a term  is  smaller  than  this  by  a factor  which 
is  the  power  of  10  equal  to  the  number  of  decimal  places  in  the  computer  word  size,  it 
cannot  affect  the  sum.  We  ignore  the  fact  that  a sum  of  small  terms  might  be  significant. 
This,  then,  defines  the  truncation  point.  Let  m be  the  number  of  the  largest  term  in  the 
sum,  k the  number  of  terms  to  be  used,  and  h the  number  of  decimal  digits  in  the  machine 
word.  Then  for  a given  k,  the  largest  absolute  value  of  the  argument  z that  can  be  used  to 
compute  g'  is  given  as 

|z3|mdm  = lz3|kdk  • 10h  . (58) 

The  power  of  ten  can  be  replaced  by  2 raised  to  a power  of  the  number  of  binary  bits  in  the 
computer  word  if  preferred.  The  coefficient  d of  eq  (57)  is  used.  Each  of  the  other  three 
should  also  be  tried,  to  find  the  smallest  number  of  the  group  for  a given  k.  Equation  (58) 
can  be  solved  for  |z|,  giving 

log  |z|  = (log  dfc  - log  dm  + h)/(3m  - 3k)  . (59) 
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A simple  computer  program  given  in  appendix  C will  find  |z|  for  each  value  of  k from  1 up 
to  the  maximum  number  of  terms  desired.  The  largest  term,  m,  is  easily  determined  because 
from  one  k to  the  next  m will  remain  the  same  or  increase  by  1 , so  it  is  only  necessary  at 
each  step  to  check  term  m + 1 to  see  if  it  is  larger  than  term  m. 

The  FORTRAN  subroutine  HANKEL  given  in  appendix  A uses  the  above  power 
series  method  to  compute  hj  and  h2  for  small  arguments.  The  coefficients  a,  b,  c,  and  d are 
given  in  lists  by  that  name.  The  truncation  points  are  given  in  the  list  called  ZMLA2.  which 
lists  values  of  |z|2  determined  by  eq  (59)  or  the  three  similar  equations. 

ASYMPTOTIC  SERIES  EXPANSION  USING  CONTINUED  FRACTIONS 

When  the  argument  z falls  outside  the  curve  in  figure  4,  h j and  h2  can  be  computed 
more  efficiently  or  more  accurately  by  asymptotic  series  than  by  power  series  methods. 
Reference  1 1 gives  information  on  branch  cuts  and  regions  of  validity  of  the  two  forms  of 
the  asymptotic  solution  (Stokes’  phenomenon).  Here  we  will  give  computing  formulas  that 
comply  with  these  requirements,  without  discussing  them  further. 

Since  a given  expansion  is  valid  in  one  or  more  quadrants,  we  choose  complete  quad- 


rants as  regions.  For  z in  quadrants  1 . 3,  or  4 use 

h2  (z)  exp  (5ir  i/12)  F2  (z)  (60) 

h2  (z)  ~ exp  (-it  i/12)G2  (z)  (61) 

For  z in  quadrant  2 use 

h2  (z)  ~ exp  (5w  i/12)  F2(z)  + exp  (1  lir  i/12)  Fj  (z)  (62) 

h2  (z)~  exp  (-ir  i/12)  G2(z)  + exp  (-7ir  i/1 2)  G|  (z)  (63) 

For  z in  quadrants  1 , 2.  or  4 use 

hj  (z)  ~ exp  (-Sir  i/12)  Fj  (z)  (64) 

hj  (z)  ~ exp  (ir  i/1 2)  G j (z)  (65) 

For  z in  quadrant  3 use 

hj  (z)  ~ exp  (- 5w  i/12)  F|  (z)  + exp  (-1  lir  i/12)  F2  (z)  (66) 

hj  (z)  - exp  (w  i/12)  Gj  (z)  + exp  (7ir  i/12)  G2  (z)  (67) 

The  four  auxiliary  functions  follow: 

M 

F,  (z)  - Kz~1/4exp(2iz3/2/3)  S CMXm  (68) 

m=0 


(69) 


M 

F2(z)  = k z’1/4exp(-2iz3/2/3)  £ CmY»> 

m=0 

M 

Gj  (z)  = k z1/4  exp  (2i  z3^2/3)  Dm  Xm  (70) 

m=0 

M 

G2  (z)  = kz^4exp(-2  iz3/2/3)  £ DmYm  (71) 

m=0 

where  X and  Y equal  * i z~3//2  respectively,  and 

K = 2,/331''6ir-1/2  = 0.853  667  218  838  951 
The  coefficients  Cm  and  Dm  are  again  computed  by  recursion  relations  where  Cq  = Dq  = 1 : 

cm  = cm-l  19  (2m-l)2  -4]/48m  (72) 

and 

Dm  = Dm_,  [9  (2m-l)2  - 16]/48m  . (73) 

Square  roots  of  z are  to  be  taken  so  that  the  real  part  of  the  root  is  always  positive  and  the 
imaginary  part  has  the  same  sign  as  the  imaginary  part  of  z.  This  applies  also  to  fourth 
roots.  The  three-halves  power  is  obtained  as  the  product  of  z and  its  square  root. 

The  summations  in  eq  (68)  -(71)  can  be  done  as  indicated  or  evaluated  by  contin- 
ued fractions.  When  done  as  indicated  they  are  asymptotic  series,  and  care  must  be  taken  to 
truncate  them  at  the  term  of  smallest  magnitude,  if  this  term  is  reached,  because  adding 
more  terms  will  reduce  the  accuracy.  Since  the  largest  term  in  these  series  will  always  be  1 , 
the  series  can  be  truncated  if  the  terms  become  less  than  10“^  in  magnitude,  where  h is  the 
number  of  decimal  digits  in  the  computer  word. 

Continued  Fraction  Expansion 

The  method  of  continued  fractions  is  more  effective  in  evaluating  these  asymptotic 
series,  and  it  is  used  in  subroutine  Hankel  in  the  FORTRAN  program  in  this  report.  The 
coefficients  are  stored  in  lists  entitled  C4,  C5,  D4,  and  D5.  In  the  remainder  of  this  section 
the  continued  fraction  technique  is  presented,  along  with  the  method  of  determining 
coefficients. 

The  continued  fraction  has  the  form 
F (x)  = b0  + aj 

x + bj  + a2  (74) 

x + b2  + T 
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It  is  to  be  used  to  evaluate  a polynomial 
M 

P(x)  - £ CmXm  . 


(75) 


m*0 


This  polynomial  can  represent  any  of  eq  (68)  -(71).  One  of  three  standard  forms  for  con- 
tinued fractions,  this  form  is  used  because  it  has  two  coefficients  at  each  stage  and  therefore 
is  equivalent  to  an  asymptotic  series  of  twice  as  many  terms.  This  reduces  by  half  the  num- 
ber of  divisions  required.  Since  complex  divisions  are  lengthy,  requiring  six  real  multiplica- 
tions and  two  divisions,  this  is  the  only  standard  form  of  the  continued  fraction  that  can 
compete  in  computer  time  with  the  asymptotic  series. 

The  coefficients  aj  and  bj  of  eq  (74)  must  be  determined  from  the  coefficients  Cm. 
The  usual  technique  is  to  express  P as  a rational  function,  then  use  the  continued  fraction  to 
evaluate  the  rational  function.  The  determination  of  the  coefficients  can  be  done  in  these 
two  steps  or  by  a second  method  which  goes  directly  from  power  series  to  continued  frac- 
tion coefficients.  The  second  method  is  preferable  because  the  loss  of  accuracy  is  more  in 
the  first.  But  since  the  first  method  is  more  easily  understood,  each  method  will  be  given;  a 
computer  program  is  included  in  appendix  C which  will  determine  coefficients  by  the  sec- 
ond method. 

Let  M in  eq  (75)  be  an  even  number  so  that  2N  = M.  (An  additional  unnecessary 
term  of  the  series  caq  always  be  used.)  The  rational  function  will  have  the  form 


N / N 

R(x)  « k ^ ®i  x/  X x*  * 
i=0  / i=0 


(76) 


where  eg  = fo  = 1 and  k = Cg.  The  coefficients  ej  and  fj  are  evaluated  from  a set  of  linear 
equations  which  can  be  described  by  displaying  a particular  case.  For  N = 3 they  are  as 
follows: 


-1 

0 

0 

Co 

0 

0 

A 

ke, 

Cl 

0 

-1 

0 

C1 

Co 

0 

A 

ke2 

c2 

0 

0 

-1 

C2 

Cl 

Co 

A 

ke3 

c3 

0 

0 

0 

C3 

C2 

Cl 

A 

fl 

c4 

0 

0 

0 

C4 

C3 

C2 

A 

f2 

c5 

_0 

0 

0 

C5 

C4 

C3_ 

A 

_f3_ 

_C6_ 

(77) 


With  ej  and  fj  thus  determined,  R(x)  is  equivalent  to  P(x)  through  the  first  M + 1 terms. 
R(x)  can  now  be  evaluated  exactly  (except  for  round-off  error)  using  a continued  fraction 
of  the  form  F(x)  of  eq  (74). 
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Rather  than  R(x),  however,  a similar  expression  in  y - l/x  is  the  form  that  is  well 
suited  to  evaluating  asymptotic  series.  This  expression  is  obtained  by  dividing  each  term  of 
R(x)  by  XN.  The  order  of  the  coefficients  can  now  be  reversed  and  a simple  algebraic  oper- 
ation can  yield  a value  of  I for  each  of  the  two  initial  coefficients  and  a new  value  for  k.  We 
will  call  this  new  rational  function  with  renamed  coefficients  ^(y).  It  will  have  the  torm  ot 
eq  (76)  but  different  coefficients,  say  e and  f instead  of  8 and  f. 

The  coefficients  aj  and  bj  are  determined  from  ej  and  fj  by  a recursive  lormula  which 
involves  constructing  an  n X n triangular  matrix  Q with  elements  q»j  as  follows: 

b0  * c0 

qj.i  * - cq  fp/aj  ‘ * 1.2,  ...,N  , 

where qj  | ■ 1.  giving  a j,  and 

bl  - fj-q|,2  • 


The  second  row: 

q2,|  * (fj-dlj+l  _bl  ql,i)/a2 


i - 2,3 N 


where  q2  2 * * • 8*vin*  a2*  and 
b2  * **1,2”  q2,3  * 

Elements  outside  the  matrix  are  assigned  a value  of  zero.  The  remaining  rows  for  m “ 3 to 


N are  as  follows: 

Mm.i  " (<lm-2,r  %i-l.i+l  ’bm-l  qm-l,i)/am 


i * m,  m+1 N , 


where  qm  m » I , giving  am,  and 

bm  * qm-l,m  ~qm.m+l  • 


The  second  method  determines  the  continued  fraction  coefficients  aj  and  bj  directly 
from  the  asymptotic  series  coefficients  Cj.  This  method  is  preferable  to  the  first  because  the 
loss  of  accuracy  in  inverting  the  matrix  in  eq  (77)  can  be  more  than  the  loss  in  this  second 

method. 


It  has  been  pointed  out*  that  the  second  method  is  probably  a variant  of  the 
Viskovakoff  algorithm  described  by  Khovanskii  (ref  12)  and  as  such  is  unstable  - subject  to 
accumulation  of  errors.  However,  it  is  sufficiently  stable  to  obtain  the  required  coefficients. 


•Privale  communication  with  AN  Stokes.  CS1RO,  Wembley . Western  Australia 

12.  The  Application  of  Continued  Fractions  and  (heir  Generalizations  to  Problems  of  Approximate  Analy- 
sis, by  AN  Khovanskii;  a monograph  in  Russian,  1956. 
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The  coefficients  are  derived  as  follows.  The  well-known  recursive  relations  that  give 
the  Nf/i  stage  of  a continued  fraction  as  a rational  function  are  used  (ref  13). 

FN(y)  = AN(y)/BN(y)  , (78) 

where 

N 

an  = X ei  yl 

i=0 

= (y  + bN)  An_,  + aN  An_2  (79) 

N 

BN  = Sfi  y1 

i=0 

= (y  + bjsj)  Bn_i  + ajsj  Bn_2  » (80) 

in  which  A_i  = 1,  Aq  = bg,  B_j  = 0,  and  Bq  = 1.  Again  y = 1/x.  The  long  division  indicated 
in  eq  (78)  is  then  carried  out,  giving  a quotient  in  terms  of  aj,  bj,  and  y that  can  be  equated, 
term  by  term,  to  the  first  2N- 1 tenns  of  the  asymptotic  series. 

The  long  division  is  carried  out  with  Ay>j  and  Bn  written  in  descending  powers  of  y. 
The  quotient  is  then  in  descending  powers  of  y or  ascending  powers  of  x.  Fortunately,  the 
first  2N+ 1 terms  determined  for  any  N are  identical  to  the  same  initial  terms  for  any  larger 
value  of  N.  This  will  be  proven  later.  The  first  few  equations  obtained  from  the  division  are 
as  follows: 

b0  = c0 

al  = C1 
-a,  b,  = C2 

al  (bl  ~a2^  = C3 

2 a2  b j -b^  + a2b2^  - C4  (81) 

From  these  equations  aj  and  bj  can  be  determined,  since  the  coefficients  Cj  are 
known.  However,  a simpler  method  is  available. 

The  long  division  indicated  in  eq  (78)  can  be  carried  to  2N+ 1 valid  places;  but  be- 
yond N+l  places,  terms  from  the  original  dividend  are  no  longer  entering  the  remainder. 
Therefore  terms  in  the  later  part  of  the  quotient  have  a simplified  form.  Since  term  n+ 1 of 


13.  Handbook  of  Mathematical  Functions  with  Formulas,  Graphs,  and  Mathematical  Tables,  ed  by  M 
Abramowitz  and  IA  Stegun;  National  Bureau  of  Standards  Applied  Mathematics  Series,  vol  55,  p 19, 
1964. 
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the  quotient  is  equal  to  the  nth  asymptotic  coefficient,  designate  it  Cm.  Note  that  the  C’s 
are  numbered  from  0 to  N.  Let  the  coefficient  of  ym  in  Bn  be  Bn  m.  Then 

N 

CN+j  " “ ^ BN.N-icN+j-i  * Kj<N  . (82) 

i*l 


Here  the  C’s  are  numerical  constants.  The  unknowns,  the  a’s  and  b’s.  are  in  the  terms  of  B. 
Suppose  that  these  unknowns  have  been  determined  up  to  n = N- 1 . Then  eq  (82)  will 
contain  two  unknowns,  ajsj  and  bjsj.  By  using  eq  (82)  for  j = N - 1 and  N,  the  unknowns  can 
be  evaluated.  The  index  N can  then  be  increased  by  I and  the  process  repeated.  The 
process  can  start  with  N * 2 if  aj , bp,  and  bj  are  provided,  but  these  are  easily  determined 
from  eq  (8 1 ).  The  terms  of  Bjq  are  determined  from  eq  (80),  which  gives  for  each  term 


Bn,m  ~ ®n-l,m-l  + Bn  Bn-J,m  + an  Bn-2,m  • 


(83) 


Any  Bn  m is  zero  if  m is  greater  than  n. 


When  j = N-l  is  used  in  eq  (82)  in  the  process  described  above,  the  coefficient  of 
the  (2N-1 ) power  of  x is  being  evaluated.  This  term  is  expected  to  contain  ajsj  aiK*  bjq.  but 
- as  will  be  proven  later  - because  the  coefficient  of  bjq  is  zero,  ajsj  is  the  only  unknown  in 
a linear  equation  and  can  be  easily  evaluated.  The  next  term  determined  with  j = N contains 
ajq  and  b^.  but  now  only  bjg  is  unknown  and  is  easily  evaluated. 


As  an  example,  theCn’s  through  n = 10  are  listed  in  table  I.  These  are  the  asymp- 
totic series  coefficients  given  by  eq  (72).  The  corresponding  an’s  and  bn’s  as  determined 
above  are  also  listed.  A more  complete  list  of  the  a’s  and  b’s  can  be  obtained  from  the 
FORTRAN  program  in  appendix  C. 


Table  1 . Asymptotic  series  coefficients,  Cn.  and  the  corresponding 
continued  fraction  coefficients,  an  and  bn. 


n 

cn 

an 

bn 

0 

1. 

0 

1. 

1 

0.10416 

0.10416 

-0.80208 

2 

0.08355 

-0.58764 

-2.28555 

3 

0.12823 

-2.29072 

-3.77864 

4 

0.29185 

-5. 11525 

-5.27462 

5 

0.88163 

-9.0628S 

-6.77193 

6 

3.32141 

7 

14.99576 

8 

78.92301 

9 

474.45154 

10 

I 

3207.49009 

A FORTRAN  program  to  compute  the  continued  fraction  coefficients  for  the  series 
given  by  eq  (72)  is  given  in  appendix  C.  This  program  can  be  easily  modified  to  determine 
the  other  set. 

Two  Proofs 

In  this  section  proofs  will  be  given  of  two  facts  used  in  the  previous  section.  Follow- 
ing this,  the  number  of  terms  required,  the  accuracy,  and  similar  topics  will  be  discussed. 

To  prove  that  the  first  2N+1  terms  of  the  quotient  Aj^/B^  are  equal  to  the  same  terms 
when  N is  a larger  integer,  use  long  division  on  eq  (79)  and  (80)  to  obtain 

An/bn  = aN-1/BN-1  + an  (AN-2  BN-1  " AN-1  BN-2)/(BN  BN-1)  • (84> 

If  the  first  quotient  on  the  right  is  to  have  terms  equal  to  the  quotient  on  the  left  up 
through  term  2N- 1 , the  remainder  must  have  no  terms  with  y to  a higher  power  than 
-(2N-1).  The  final  divisor,  Bfsj  Bn_j,  contains  y to  the  (2N-1)  and  lower  powers.  There- 
fore, the  proof  is  complete  if  the  numerator  of  the  remainder  is  a constant.  To  show  this, 
use  eq  (79)  and  (80)  to  evaluate  Bjq-1  and  A^j-i ; it  can  be  shown  that 

aN-2  BN-1  " AN-1  BN-2  = _aN-l  (AN-3  BN-2  " AN-2  BN-3) 

= (-1)N  aN-l  aN-2  • • • al  • 

The  right-hand  product  is  obtained  by  repeatedly  applying  the  middle  result.  The  product 
of  a’s  is  a constant,  completing  the  proof. 

The  second  proof  required  is  that  in  the  quotient  of  Ajq/Bj^j,  C2N  (the  coefficient  of 
y-2N)  wiii  contain  no  aj  or  bj  to  higher  than  term  N and  C2N+1  (the  coefficient  of  y~2N-l ) 
will  involve  no  aj  to  higher  than  term  N+ 1 and  no  bj  to  higher  than  term  N. 

The  first  part  is  intuitively  obvious.  Since  from  the  preceding  proof  C2N  will  be  the 
same  when  derived  from  the  ratio  Ax/Bx  for  any  x as  long  as  it  is  N or  greater,  we  need 
consider  only  the  case  where  x is  N.  But  since  from  eq  (79)  and  (80)  Ajsj  and  Bj»j  contain  no 
a’s  or  b’s  of  greater  than  term  N,  C2N  cannot  contain  a’s  or  b’s  of  higher  terms. 

By  the  same  argument  C2N+1  can  contain  no  a’s  or  b’s  to  higher  terms  than  N + 1. 
There  remains  to  be  proven  only  that  bjq+j  cannot  exist  in  C2N+I  or  that  its  coefficient, 
which  we  will  call  E,  is  zero.  Applying  eq  (82)  for  N+ 1 and  j = N gives 

N+l 

CN+1  + N = ‘ X BN+l,N+l-i cN+l+N-i  • 
i=l 

E,  the  coefficient  of  bn+  j in  this  expression  from  eq  (83),  takes  the  form 
N+l 

BN,N+l-i  C2N+l-i  • 
i=l 
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But  by  choosing  j = N in  eq  (82)  we  see  that  terms  1 to  N for  C 2N  equal  terms  2 to  N+ 1 in 
E,  so 

E = ~BN,NC2N  + C2N  • 

However,  since  Bjq  jq>  the  coefficient  of  yn  in  Bn,  is  always  1 by  eq  (80),  E = 0.  Therefore 
bjq+j  does  not  exist  in  Ctn+j. 

Number  of  Terms 

The  number  of  terms  or  stages  to  use  in  the  continued  fraction  was  arrived  at  by  a 
trial  and  error  process.  For  a given  number  of  terms,  a real  positive  argument  was  decreased 
until  the  accuracy  began  to  drop.  The  magnitude  just  before  this  drop  was  considered  to  be 
the  optimum  point  to  increase  the  number  of  stages  by  one.  Because  the  argument  to  the 
continued  fractions  is  z^/2,  we  took  the  larger  of  the  magnitudes  of  the  real  and  imaginary 
parts  of  z^/2  as  the  test  number.  This  number  is  then  compared  to  the  3/2  power  of  the 
points  determined  along  the  real  axis  by  trial  and  error. 

The  above  method  appears  to  work  well  although  it  involves  no  thorough  under- 
standing of  the  way  complex  numbers  affect  the  successive  convergents  of  a continued 
fraction.  Table  2 shows  the  points  down  to  which  a given  number  of  stages  gives  full  accu- 
racy for  positive  real  arguments  and  lists  the  3/2  power  of  these  numbers  as  used  in  the 
FORTRAN  program  list  called  ZMLA5. 

Division  Lines 

The  power  series  method  is  now  to  be  used  for  small  arguments  and  the  continued 
fraction  method  for  large  arguments.  The  exact  dividing  line  between  them  is  needed.  The 
division  line  of  figure  4 was  arrived  at  by  computing  the  functions  along  rays  from  the  ori- 
gin, using  both  power  series  and  continued  fractions.  The  number  of  decimal  places  to 
which  the  functions  determined  by  the  two  methods  agree  tends  to  reach  a maximum  at 
some  distance  from  the  origin  along  each  ray.  At  distances  short  of  this  maximum  we  can 
assume  that  the  continued  fraction  method  is  less  accurate  than  the  power  series.  At  dis- 
tances beyond  the  maximum,  the  power  series  is  assumed  to  be  less  accurate.  The  maxi- 
mum therefore  indicates  the  ideal  place  to  change  from  one  method  to  the  other  if  the 
decision  is  to  be  based  solely  on  accuracy.  This  method  was  used  to  determine  figure  4. 

A complication  arises,  however.  Along  certain  rays  from  the  origin,  b]  and  its  deriva- 
tive reach  a maximum  number  of  places  at  very  different  distances  from  h2  and  its  deriva- 
tive. The  principal  problem  is  at  ±60°  but  persists  from  about  30°  to  90°.  At  60°,  hj  is 
small  in  magnitude  and  h2  is  large.  The  power  series  method  cannot  compute  the  small 
values  accurately  due  to  loss  in  accuracy  in  subtraction  in  eq  (48).  The  accuracy  of  the 
continued  fraction  for  h2  is  poor  at  60°  because  eq  (69)  becomes  a nonalternating  series  and 
continued  fraction  approximations  are  not  known  to  improve  the  accuracy  of  nonalternat- 
ing asymptotic  series  as  they  do  for  alternating  series. 

A reasonable  solution  to  this  problem  is  to  compute  hj  by  continued  fractions  and 
h2  by  power  series  for  arguments  at  these  angles  and  magnitudes  from  4 to  10.  However,  as 
will  be  shown  later,  the  above  solution  has  not  been  employed  at  this  time  since  this  area  is 
not  of  great  importance  for  normal  mode  computations.  Instead,  the  argument  was  chosen 
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Table  2.  Cut-off  points  for  determining  the  number  of  stages  in  the  continued  fractions. 


Number  of  Stages 

Real  Argument 

X 

Progrra  Test  Value 
x3/2 

! 

106 

109 

2 

80 

715.0 

3 

35 

207.0 

4 

22 

103.0 

5 

13 

47.0 

6 

11 

36.4 

7 

9 

27.0 

8 

8 

22.6 

9 

7 

18.5 

10 

6.5 

16.6 

11 

6 

14.7 

12 

5.8 

14.0 

13 

S.S 

12.9 

14 

5.3 

12.2 

IS 

5.1 

11.5 

16 

4.9 

10.8 

17 

4.5 

9.5 

18 

4.4 

9.2 

that  gave  equivalent  accuracy  for  the  two  methods.  Along  60°  this  minimum  accuracy  is  9 
decimal  places. 

The  following  relationship  exists  between  hj  and  h2  for  positive  and  negative  values 
of  the  imaginary  part  of  the  arguments: 

h,  (z*)  - lh2(z»*  , 

where  the  * means  complex  conjugate.  Thus,  the  above  discussion  at  60°  can  be  translated 
to  -60°.  Also,  the  functions  actually  need  to  be  computed  only  in  quadrants  I and  II.  They 
could  then  be  evaluated  in  quadrants  III  and  IV  by  the  above  relationship.  The  above  rela- 
tionship explains  the  symmetry  of  figure  4 about  the  real  axis. 
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COMPARATIVE  ACCURACY 


The  accuracy  of  the  three  methods  - power  series,  asymptotic  series  and  continued 
fractions  — has  been  determined  on  a C DC  computer  with  48  bits  or  14.4  decimal  places  of 
accuracy  in  the  floating  point  word.  Since  this  differs  from  the  double  precision  word 
length  of  60  bits  or  1 8. 1 decimal  places  that  applies  to  the  preceding  part  of  this  report, 
these  results  arc  for  comparative  and  illustrative  purposes  only. 

Accurac'  is  determined  by  computing  the  functions  and  either  comparing  the  an- 
swers for  the  se’  eral  different  computing  methods  or  computing  the  wronskian.  The 
wronskian  is  a constant  given  by  the  relationship 

hj  h'2-h2h',  = - 1.45 7495441 04 i = -i96J/3/7r  . (85) 

The  wronskian  will  determine  the  accuracy  of  the  functions  if  it  can  be  computed  without 
loss  of  accuracy.  If  the  two  products  in  it  are  large,  though,  accuracy  will  be  lost  in  the 
subtraction.  This  generally  happens  for  arguments  near  the  negative  real  axis.  Here  accura- 
cy must  be  determined  by  comparing  answers  from  different  methods.  The  accuracy  of  the 
functions  and  their  derivatives  will  generally  be  about  equal. 

Figure  5 illustrates  the  accuracy  that  is  obtained  in  different  parts  of  the  complex 
plane  of  the  argument,  z,  by  using  the  power  series  method.  On  the  inner  contour,  the 
functions  hj  and  h2  and  their  derivatives  have  12  places  of  accuracy.  On  the  outer  contour, 
the  accuracy  is  1 1 places.  As  expected,  the  accuracy  is  best  for  arguments  of  small  magni- 
tude. The  accuracy  remains  best  in  directions  from  the  origin  in  which  the  functions  are 
large  in  magnitude.  This  is  because  less  accuracy  is  lost  in  subtraction.  Accuracy  must  be 
lost  when  individual  terms  of  the  series  are  large  but  the  sum  is  small. 

Figure  6 shows  accuracy  contours  for  the  asymptotic  expansion  with  both  the  direct 
and  continued  fraction  evaluation  of  the  series.  Here,  the  best  accuracy  is  obtained  for  large 
arguments,  and  accuracy  decreases  toward  the  origin.  As  can  be  seen,  each  of  the  two  meth- 
ods is  better  in  some  directions  from  the  origin.  The  choice  of  methods  then  depends  upon 
which  directions  are  of  most  value  to  the  normal  mode  program.  The  dots  on  the  figure 
show  the  locations  at  which  the  functions  were  evaluated  in  a typical  surface  duct  run. 
Although  arguments  can  lie  anywhere  in  the  plane,  most  of  them  follow  this  pattern.  They 
lie  just  above  the  negative  real  axis  and  in  a narrow  angle  above  the  positive  real  axis.  The 
continued  fraction  method  is  distinctly  better  on  this  positive  side.  Since  computing  time 
also  favors  the  continued  fraction  method,  it  is  clearly  the  method  to  use. 

If  the  12-place  accuracy  contour  from  figure  6 lies  inside  that  for  figure  5 at  some 
angle  from  the  origin,  12  places  can  be  obtained  at  any  range  along  this  angle  by  using  either 
power  series  or  asymptotic  expansion  in  the  interval  of  overlap.  If  the  asymptotic  expan- 
sion contour  lies  outside  the  other,  there  is  an  interval  in  which  1 2 places  cannot  be  ob- 
tained. Only  some  lesser  number  of  places  can  be  obtained  in  this  interval.  These  contours 
apply  when  both  functions  and  their  derivatives  are  all  computed  by  a single  method.  As 
mentioned  earlier,  increased  accuracy  could  be  obtained  in  some  areas  by  computing  the 
two  functions  by  different  methods. 
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Figure  S.  Locus  of  arguments  for  which  the  power  series  evaluation  of  the 
modified  Hankel  functions  gives  12  and  1 1 decimal  places  of  accuracy  for 
a computer  word  length  of  14.4  decimal  places. 


Figure  6.  Locus  of  arguments  for  which  the  direct  and  continued  fraction 
evaluation  of  the  asymptotic  series  gives  1 1 and  1 2 decimal  place  accuracy. 
The  arguments  at  which  the  modified  Hankel  functions  were  evaluated  in  a 
typical  normal-mode  run  are  shown. 
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MODE  FOLLOWER  PROGRAM 


Appendix  D lists  the  Mode  Follower  Program  in  FORTRAN.  It  is  not  a part  of  the 
general  normal  mode  program,  but  is  related  in  that  it  uses  some  parts  of  the  general  pro- 
gram. The  purpose  of  the  mode  follower  is  to  trace  a given  eigenvalue  as  some  parameter  is 
varied.  This  parameter  is  usually  frequency,  but  any  profile  parameter  can  also  be  varied. 

The  eigenvalues  at  a given  set  of  parameters  are  discrete  points.  By  permitting  the  parame- 
ter to  vary,  the  eigenvalues  become  a set  of  lines,  and  this  often  clarifies  their  behavior  at  the 
fixed  points.  Figures  7-9  illustrate  this. 

Figure  7 is  a sound  speed  profile  consisting  of  two  ducts.  Figures  8 and  9 show  the 
real  and  imaginary  parts  of  some  eigenvalues  of  the  profile  over  a range  of  frequencies.  The 
imaginary  parts  are  expressed  as  mode  attenuations.  The  figures  show  a region  where  both 
ducts  are  exerting  an  influence  on  the  eigenvalues.  The  broken  lines  show  the  location  of 
eigenvalues  for  a profile  that  consists  of  only  the  upper  duct  of  figure  7.  Considerable  time 
could  be  spent  studying  the  interaction  between  the  two  ducts,  but  since  the  purpose  here  is 
to  illustrate  eigenvalues  as  functions  of  a parameter,  only  a brief  description  will  be  given. 

Modes  are  numbered  by  the  real  parts  of  their  eigenvalues.  This  numbering  is  con- 
sistent with  the  number  of  beats  or  changes  of  n in  the  phase  of  the  depth  functions.  Thus 
the  eigenvalue  of  a mode  numbered  1 in  a profile  consisting  of  only  the  upper  duct  lies 
exactly  over  the  eigenvalues  of  a mode  in  the  double  duct  in  figures  8 and  9,  but  this  mode 
in  the  double  duct  changes  number  each  time  it  crosses  the  real  part  of  another  mode.  The 
depth  function  actually  gains  an  additional  beat  each  time  this  happens.  The  background  of 
modes  that  are  being  crossed  consists  of  the  higher  order,  untrapped  modes  associated  with  the 
lower  duct. 

Mode  2,  of  the  upper  duct  only,  does  not  have  a single  mode  in  the  double  duct  that 
overlies  it  exactly.  Instead,  a mode  attempts  to  follow  it  at  frequencies  above  1350  Hz. 

Below  this  frequency,  successive  modes  follow  its  path  for  short  intervals.  This  interplay 
between  modes  occurs  when  mode  2 of  the  upper  duct  is  in  some  sense  equally  as  untrapped 
as  the  modes  associated  with  the  lower  duct. 

The  imaginary  parts  of  the  modes  follow  similar  patterns;  but  because  the  mode 
numbering  is  not  determined  by  the  imaginary  parts,  the  mode  numbers  sometime  jump 
from  one  line  to  another.  An  important  feature  of  these  two  plots  is  that  if  the  real  parts  of 
the  eigenvalues  cross,  the  imaginary  parts  do  not;  and  vice  versa.  Thus  two  eigenvalues  do 
not  tend  to  become  equal  at  a point  which  would  make  them  degenerate. 

The  mode  follower  program  will  tend  to  follow  the  continuous  curves.  Thus  if  start- 
ed in  the  right  direction  on  mode  59  at  1450  Hz,  it  will  follow  along  the  continuous  mode 
which  becomes  successively  mode  58,  57,  56,  and  55. 

Figures  such  as  8 and  9 can  be  drawn  by  computing  the  eigenvalues  at  a sufficient 
number  of  frequencies  to  determine  the  lines.  The  mode  follower  does  this  for  a given 
eigenvalue  while  adjusting  the  step  size  so  the  mode  will  not  be  lost,  or  so  the  program  will 
correctly  follow  the  mode.  The  step  size  is  permitted  to  become  large  where  the  eigenvalue 
can  be  approximated  by  a parabolic  curve,  but  it  shortens  when  extrapolation  to  the  next 
point  becomes  less  accurate. 
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Figure  9.  The  imaginary  parts  of  the  modes  whose 
real  parts  are  shown  in  figure  8 expressed  as  mode 
attenuation.  To  avoid  confusion,  they  are  not 
shown  across  the  full  frequency  interval. 

When  frequency  is  the  variable  parameter,  the  group  velocity  of  the  mode  can  be 
computed  easily  since  a numerical  derivative  can  be  computed.  Group  velocity  is  given  by 
the  relationship 

Cg  = dto/d  (Re  k) 

3:  Acu/A  (Re  k) 

3s  - A f v^/f  Av  , 

where  k is  the  horizontal  wave  number  of  the  mode  and  v is  the  real  phase  velocity.  The 
mode  follower  prints  this  value  out  at  each  step,  along  with  the  eigenvalue. 

IMPLEMENTATION  OF  THE  MODE  FOLLOWER 

The  mode  follower  was  originally  implemented  for  a two-layer  normal  mode  which 
differed  from  the  n-layer  program  in  that  the  derivative  dG/dv  of  the  characteristic  equation 
was  evaluated  along  with  G.  The  iteration  for  roots  of  G was  thus  Newton-Raphson  and  is 
given  by  the  relationship 


vi+l  * Vj-G/G'  . (86) 

This  is  simpler  than  the  secant  iteration  of  eq  (15),  in  which  G'  must  be  evaluated  numeri- 
cally. Because  of  the  simpler  iteration,  an  effective  scheme  for  mode  following  was  available. 


Since  considerable  effort  was  required  to  develop  a similar  scheme  for  the  n-layer  case,  the 
two-layer  mode  follower  will  be  briefly  described  to  serve  as  an  introduction  to  the  n-layer 


The  two-layer  mode  follower  employs  one  iterative  step  of  eq  (86)  at  each  point 
where  G is  evaluated.  Thus,  a root  that  is  inexact  but  sufficiently  exact  is  obtained.  The 
original  estimate  is  obtained  by  extrapolating  from  the  three  most  recent  roots.  If  this  esti- 
mate is  sufficiently  close  to  the  true  root,  the  single  iterative  step  will  make  a small  correc- 
tion. G/G',  that  will  bring  the  estimate  very  close  to  the  true  root.  By  using  the  size  of  this 
correction  to  control  the  step  size,  the  program  is  self-regulating.  The  program  works  well 
when  a permissible  value  of  G/G'  of  10-6  to  I0~4  is  used.  Outside  this  interval  the  step  size 
is  either  doubled  or  halved. 

The  multiple-layer  program  differs  from  this  in  several  details.  The  extrapolation 
from  the  previous  three  points  is  done  not  only  for  the  phase  velocity  but  also  for  the  nu- 
merical derivative, 

D*1  = Av/AG. 

Lagrange  three-point  interpolation  is  used,  given  by  the  form 

v(x  j)  (x  - x2>  (x  - X3)  v(X2>(x  - X|)(x- X3)  v(X3)(x-xj)(x-X2> 

V(X)  * (X(  — x2)(xj  -x3)  ' (x,  - x2)(X2  - X3)  + (X|  - X3MX2  - X3)  * 

(87) 

where  x is  the  new  value  of  the  parameter  that  is  being  varied  (usually  frequency)  and  X| . 

X2.  and  X3  arc  the  three  previous  values,  xj  being  the  most  recent.  To  extrapolate  the 
derivative,  v is  replaced  by  D-'  in  eq  (87).  Both  quantities  are  complex  numbers. 

The  determinant  is  now  evaluated  at  this  new  phase  velocity  to  give  a value  Gq 
Next  a corrected  value  of  phase  velocity,  vq.  is  obtained: 

v0  M v-Gq  D-1  . (88) 

In  the  two-layer  case,  the  size  of  the  correction.  GO-* , is  used  (0  control  the  step  size. 
Because  the  numerical  derivative  is  less  precise,  we  evaluate  G once  more  at  this  new  posi- 
tion. obtaining  G | . A new  numerical  derivative  is  next  calculated : 

D0~ 1 * (vo  “ vV(Gi  - Gq)  . 

This  derivative  is  now  compared  with  the  extrapolated  value  to  determine  whether  the  step 
size  should  be  changed.  To  do  this  an  error 

E - II  -D0/D|2 

is  computed.  Good  results  have  been  obtained  by  keeping  E between  10“^  and  10“^.  If  E 
becomes  larger  than  this,  the  step  size  is  halved  and  the  extrapolation  is  tried  again.  Should 
halving  the  step  size  five  times  fail  to  obtain  a value  of  E less  than  10“'.  the  mode  is  pre- 
sumed to  be  lost  and  the  program  halts. 
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If  E is  less  than  10“2,  the  step  is  successful  and  the  stored  values  are  updated  for  the 
next  step.  Before  v is  stored,  though,  the  iterative  step  of  eq  (88)  is  applied  one  more  time 
to  obtain  a more  precise  value  of  v.  This  requires  little  extra  effort  because  the  quantities 
Gj  and  Dq~*  are  already  available. 

If  the  error  E is  less  than  10“^,  the  next  step  size  is  doubled. 

It  is  possible  for  the  extrapolation  to  be  too  successful.  That  is,  if  v is  very  near  the 
true  root.  Gq  and  Av  will  be  very  small.  The  numerical  derivative  may  then  be  inaccurate. 
Therefore,  before  the  error  term  E is  computed,  a quantity 

F = |v/Av|2 

is  computed.  If  F is  greater  than  10*^,  the  extrapolated  derivative  is  used  rather  than  the 
computed  derivative  and  the  program  proceeds  to  the  next  step.  If  F is  greater  than  lO^4, 
the  step  size  is  doubled  before  proceeding  to  the  next  step. 

The  other  principal  part  of  the  program  is  the  initialization  which  must  evaluate  v at 
three  values  of  x to  obtain  the  numbers  needed  for  the  first  extrapolation,  eq  (87). 

INPUT  AND  OUTPUT 

The  first  input  card  contains  the  maximum  number  of  steps  allowed,  the  limits 
applied  to  E and  F,  and  keys  which  control  both  the  amount  of  detail  in  the  printout  and 
whether  the  profile  parameters  are  to  be  read  in  or  retained  from  the  previous  run.  Default 
values  are  supplied  when  these  items  are  left  blank.  Next  the  profile  parameters  are  read  in. 
These  are  an  older  style  and  only  permit  specification  of  the  absorption  loss  at  the  top  of  a 
layer.  The  sound  speed  gradient  is  assumed  to  be  real  at  the  top  of  any  layer. 

A final  card  indicates  which  variable  — frequency,  sound  speed,  depth,  gradient, 
absorption,  or  density  - will  be  varied,  by  specifying  a number  called  nx  in  the  program, 
from  1 to  6.  The  next  number,  ny,  specifies  which  layer  the  variable  will  be  in.  This  layer 
number  is  not  needed  if  frequency  is  selected.  A third  number,  nz,  indicates,  if  zero,  that 
the  profile  will  remain  continuous  as  the  selected  parameter  is  varied.  If  nz  is  not  zero,  the 
selected  parameter  moves  alone  without  a compensating  motion  in  other  profile  parameters. 
The  card  next  gives  the  initial  and  final  value  of  the  parameter  to  be  varied  and  the  initial 
step  size.  Finally,  the  particular  mode  to  be  followed  is  indicated  by  giving  an  approximate 
phase  velocity  and  an  initial  step  size.  These  must  be  chosen  such  that  the  subsequent  itera- 
tion will  converge  on  the  correct  mode. 

The  principal  output  of  this  program  is  the  print  statement  at  line  314.  Each  line  of 
output  contains  the  value  of  the  parameter  being  varied,  the  complex  phase  velocity,  the 
determinant  G,  the  derivative  D-* , the  error  term  E,  the  mode  attenuation,  the  mode  group 
velocity  (if  frequency  is  the  parameter  being  varied),  and  the  step  number.  After  the  final 
step,  the  profile  in  its  final  form  is  printed  out. 


CONCLUSIONS 


1 . An  effective  program  for  computing  propagation  loss  in  a layered  ocean  by  nor- 
mal modes  has  been  developed.  Complete  documentation  for  the  program  is  contained 
herein. 

2.  Sediment  layers  are  modeled  as  fluids  in  which  densities,  sound  speeds,  and  ab- 
sorption can  be  specified.  This  permits  a complete  wave  solution  for  bottom  reflected 
sound  energy. 

3.  A continued  fraction  technique  for  evaluating  asymptotic  series  is  shown  to  give 
superior  results  in  evaluating  the  auxiliary  functions  required  in  this  program,  the  modified 
Hankel  functions  of  order  1/3. 

4.  A mode  follower  program  given  here  is  useful  in  tracing  eigenvalues.  Such  traces 
are  needed  to  understand  the  eigenvalue  structure. 


RECOMMENDATIONS 

1 . Improve  the  mode  locating  ability  of  this  normal-mode  program  to  make  it  self- 
contained.  It  currently  requires  user  interaction  to  locate  eigenvalues. 

2.  Investigate  methods  to  incorporate  the  effect  of  rough  boundaries  into  this 
program. 
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APPENDIX  A:  NORMAL  MODE  PROGRAM  IN  FORTRAN 

This  program  consists  of  the  main  program  and  seven  subroutines.  The  main  pro- 
gram handles  the  input  and  output  and  performs  much  of  the  computation.  This  includes 
profile  preparation,  mode  search,  determination  of  depth  function  coefficients,  normaliza- 
tion, computation  of  depth  functions,  and  summation  of  modes.  Auxiliary  functions  are 
performed  by  the  subroutines  SETUP  and  DET,  which  set  up  the  determinant,  then  evaluate 
it.  This  is  the  determinant  from  which  eigenvalues  are  determined.  The  subroutine  HZERO 
determines  the  Hankel  functions  of  order  zero,  second  type,  which  gives  the  range  depend- 
ence of  the  modes.  Only  a single  term  of  the  asymptotic  expansion  is  needed  for  this 
function. 

Subroutine  HANKEL  evaluates  the  modified  Hankel  functions  of  order  1/3,  by 
which  the  depth  dependence  of  the  modes  is  expressed.  The  majority  of  computing  time  is 
usually  expended  in  this  subroutine.  Subroutine  CFR  is  used  by  subroutine  HANKEL  to 
evaluate  continued  fractions.  Subroutine  RCOEF  evaluates  and  prints  reflection  coeffi- 
cients when  they  are  requested. 


(HAUM^OtfO^lOtUIAUM^OlOQDSOlUlAiwMxOtOGD^OMII^UM^OtOGB^aUIAUM^OtfXB^aiUlAUM 
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1 C THIS  IS  THE  MAIN  PART  OF  NLAYNM 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  LAMBOA, LAMBDI 
INTEGER  COL 

REAL  R ATTEN,  T RE , RX 
DIMENSION  LOSPCH( 5 , 25) 

COMMON  /HAN/  H2R,H2I,H1R,H1I,H2PR,H2PI,H1PRfH1PI,R 
COMMON/INPUT/  Z(12),  N,  OMEGA,  V,  VI,  CON(12),  GSQ<12), 

1 CAY (12),  LAMBDA,  LAMBDI,  G(12) 

2 , RHO( 12),  GI (12),  G SO  1(12),  CAYI(12) 

COMMON  /EXPO/  EXSUM,  CNTR , RATIO(25) 

COMMON/DETMNT / A(25,4),  0(25,4) 

COMMON/PARTS/  ZT ( I 2 ) , ZT I ( 1 2 ) , ZB( 1 2 ) , ZBI ( 1 2 ) 

COMMON/ RE FL/  AF(12,200),  AG(12,200),  BF(12,200),  BG(12,200), 

2 EIGEN ( 350 ) , EIGEN  I ( 350 ) , B (25,4),  BI(25,4),  CB(12),  CBI(12). 

3 CAYSQM2).  CAYSQ [(12),  NN 

DIMENSION  D( 350 ) , DI(350),  F(100),  FI(IOO),  HZER02(350), 

* DA( 350 ) , SRES( 350 ) , GAMMA  1(12),  BLPK(12), 

4 HZER2I ( 350 ) , DPK(12),  GCU(12),  GCUI(12),  C I ( 1 2 ) . 

3 PHASE  V ( 350 ) , PHASI  V(350),  UU(2000),  UUI(2000> 

COMMON  /LIMIT/  TLIM.  EXPONT,  SLIM 
DIMENSION  LOSS( 101 ) 

DIMENSION  C( 1 2 ) , 0EPTH(52 ) , DBLOSS(350),  C0L(120), 

ICONTR(IO),  EF(2).  FMAG( 350 ) , FANG(lOO), 

2GAMMA( 12),JSMBL(10), JCOUNT ( 5 ) , JCOU( 5 ) , LEVE  L ( 41 ) ,PLEV(5 ) , RLOSS( 100) 
3 , RLOS( 101 ) , RECVRS(51 ) .TEST (3) , ING( 11) 

EQUIVALENCE  ( FF , E F ( 1 ) ) ,(DEPTH( 1 ) , SOURCE ) , ( DEPTH( 2 ) , RECVRS( 1 ) ) , 

1 ( R LOS ( 2 ) , RLOSS ( 1 ) ) 

COMMON  /AHZERO/  HZ EROR .HZEROI 

DATA  ( CONTR(I),  1=1,4)  / 1 1 0 . DO , 95 . DO , 80 . 00 , -1 000 . DO/ , 

1 (JSMBL(I),  1 = 1 , 3 ) / 1 HI , 1 H* , 1 H8/ , 

* ( ING( I ) , I = 1 ,10)/1H0,1H1 ,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/, 

2 ( TEST ( I ) , 1=1,3)  /.2D0, 1 .DO, 5. DO/ 

TLIM  = 25. 

SLIM  = -8.54988 

C READ  IN  PARAMETERS 

I READ  11,  HI , K2 , K3,  K4,  K5,  K6 , K7,  KB,  K9 
C KEYS*  1-DEPT  FN  PRINTOUT,  2-LOSS  PRINTOUT,  3-REFLECTION  COEFF  PRINTOUT 
C 4-CHANGE  CONTOURS,  5-CONTINUE  MODES 

II  FORMAT  (1014) 

PRINT  u,  M,  K2 , K3 , K4,  K5,  K6,  K7,  KB,  K9 
13  FORMAT  (6H1KEYS  , 1014) 

IF  (K1  . LT.  7)  GO  TO  5 
READ  1 1 , M 

READ  434,  (SRES(I),  1=1 ,M) 

434  FORMAT  (5D16.7) 

5  IF  (K8  .NE.  1 ) GO  TO  8 
REAO  20,  T LIM 
SLIM  • -DCBRT(TLIM**2) 

PRINT  30,  TLIM,  SLIM 
8 EXPONT  = DEXP ( (TLIM  ♦ TLIM)  / 3.) 

K6  ■ K6  + 1 
INSUR  > 0 

IF  (K2  . LT.  10)  GO  TO  16 
K2  - K2  - 10 
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16 

MPCH  • 0 

50 

1 F (45  .LT.  10)  GO  TO  17 

59 

K5  ■ 45  - 10 

60 

MPCH  ■ 1 

61 

17 

IF  (44  .NE.  1 ) GO  TO  3 

62 

READ  20.  (CONTR(I),  I ■ 1 ,9) 

63 

CONTR(IO)  » -1000.00 

64 

READ  4,  (J  SM6KI),  I *1,9) 

65 

4 

FORMAT  ( 9A1 ) 

66 

3 

READ  10.  N,  FREQ 

67 

10 

FORMAT  (12.  010.1) 

68 

IF  (N.EQ.O)  GO  TO  999 

69 

2 

PRINT  12,  N,  FREQ 

70 

12 

FORMAT  (13,  8H  LAVERS,  .F10.1,  3H  HZ) 

71 

READ  20,  (2(1).  1-1 ,N) 

72 

PRINT30 , (Z(I),  1-1, N) 

73 

READ  20,  (C(I),  I - 1.N) 

74 

20 

FORMAT  (8D10.4) 

75 

PRINT30 , (C(I),  I - 1 ,N) 

76 

READ  20,  (CB( I),  I « 1,N) 

77 

PRINT30,  ( CB( I ) , I • 1,N) 

78 

RF»n  ?f>,  (C.amMA(T),  I - 1.N) 

79 

PRINT30,  ( GAMMA ( I ) , I * 1 . N ) 

80 

READ  20,  (DPK(I),  I - 1,N) 

81 

PRINT  30,  ( DP4 ( I ) , I - 1 , N ) 

82 

READ  20.  (BLPK(I),  I - 1.N) 

83 

PR1NT30,  (BLPK(I),  I - 1 , N) 

84 

READ  20,  (RHO(I),  I - 1 . N ) 

95 

PRINT  30,  (RHO(I),  I - 1 , N) 

86 

IF  (FRED  .GT.  0.)  GO  TO  18 

87 

FREQ  - - FREQ 

88 

ATTEN  - 0. 

89 

GO  TO  19 

90 

18 

F SQ  • (FREQ  / 1000. )-*2 

91 

ATTEN  - .1  • F SQ  / (1.  ♦ F SQ)  ♦ 40. 

92 

19 

ATTEN  • ATTEN  * 1.0936 

93 

PRINT  14, ATTEN 

94 

14 

FORMAT  (8H  ATTEN  > .G10.5,  5HDB/KM  ) 

95 

ATTEN  « ATTEN  / 1000.00 

96 

30 

FORMAT  (9F14.5) 

97 

C COMPLETE  PROFILE 

98 

DO  33  I - 1 ,N 

99 

IF  (RHO(I)  .EQ.  0.)  RHO(l)  > 1.02 

100 

IF  ( CB( 1 ) .NE.  0. ) GO  TO  31 

101 

CB( I ) • C(I+1 ) 

102 

31 

IF  (C(I ) .NE.  0. ) GO  TO  32 

103 

C(I)  - CB(I-I) 

104 

32 

IF  (DPK(I)  .GE.  0. ) GO  TO  34 

105 

OKI)  • CBI(I-I) 

106 

GO  TO  36 

107 

34 

Cl ( I ) • 0. 

108 

IF  (DPK( I)  .EQ.  0. ) GO  TO  36 

109 

T * 27287.52708  / DPK(I) 

110 

CI(I)  « T - SQRT ( ( T - C( I ) ) • (T  ♦ C( 

111 

36 

CBI(I)  - 0. 

112 

IF  ( GAMMA ( I ) .NE.  0.)  GO  TO  38 

113 

C 

••BOTH  SOUND  SPEEDS  GIVEN 

* f so  / (4ioo.  ♦ f so) 
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OLDFN 
DO  41 


IF  (BLPK(I)  .IE.  0.)  GO  TO  37 
T - iuoi  . 9^/0b  / BLPK(I) 

CBI(I)  - T - SORT  ( ( T - C B ( I ) ) • (T  4 CB( I ) ) ) 

T - C(I)  • ( C ( I ) • • 2 - 3.  • Cl ( I )**2) 

T I ■ CI(I)  • (3.  • C(I)**2  - Cl ( I ) **2) 

IF  (BLPK(I)  . LT . 0. ) GO  TO  39 
TEMP  » CB( I ) * *2  - CBI ( I )**2 
TEMPI  • 2.  • CB( I ) • CBI ( I ) 

OENOM  * TEMP**2  4 TEMPI *• 2 
TEMP  « TEMP  / DENOM 
TEMPI  * -TEMPI  / OENOM 

GAMMA(l)  • 0.5  * (C(I)  - (T  • TEMP  - TI  * TEMPI))  / 

• (2(1*1)  - 2(1)) 

GAMMA  1(1)  . 0.5  • (CI(I)  - (T  « TEMPI  4 TI  • TEMP))  / 

• (2(1*1  ) - 2(D) 

IF  (I  . EQ.  N)  GO  TO  27 
GO  TO  33 

**SPECIAL  CASE,  GRADIENT  REAL  NUMBER 

IF  (Cl ( I ) . EQ . 0.)  GO  TO  42 

TEMP  « CB( I ) • *2 

TEM  « TEMP*»2 

TEMPI  . CI(I) 

C0EF1  > C I ( I > 

C0EF2  * 2.  • TEMP  • CI(I)  4 TI 
C0EF3  » 2.  * T * CB( I ) 

C0EF4  » TEM  * CI(I)  - TEMP  • TI 
OLDFN  « 1 .020 
DO  41  J > 1,10 

FN.  ( ( ( COEF 1 * TEMPI)  4 C0EF2)  • TEMPI  4 C0EF3)  * TEMPI  ♦ C0EF4 

FP  • ((4,  • COEF 1 • TEMPI)  42.*  C0EF2)  * TEMPI  4 C0EF3 

TEMPI  * TEMPI  - FN  / FP 

IF  (FN  .GE.  OLDFN)  GO  TO  43 

OLDFN  « FN 

CONTINUE 

CBI ( I ) * TEMPI 

GAMMA ( I ) • .5  »(.5  • (CI(I)  * ( TEMP  - CBI(I)**2)  - TI)  / 

* (CBI I ) • TEMPI)  4 C( I ) ) / (2(141)  - 2(1)) 

GO  TO  28 

GAMMA(I)  - C ( I ) - T / (CB( I ) **2  * (2(141)  - 2(1))  • 2.) 

GO  TO  33 

• •SO'JNP  S»FED  AND  GRADIENT  GIVEN 
IF  (I  .EQ.  N)  GO  TO  33 
T • C(I)  • (C( I )*»2  - 3.  * C I ( I ) **2 ) 

TI  ■ CHI)  • (3.  • C(  I ) **2  - CI(I)**2) 

IF  (BLPK(I)  .EQ.  0.)  GO  TO  29 
IF  (BLPK(I)  .LT.  0.)  GO  TO  2B 

TEMP  ' (BLPK(I)  / 54575. 05416)**2  / (2(141)  - 2(1))  • 0.5 
T • T • TEMP 
TI  » TI  • TEMP 

TEMP  . .5  • C(I)  / (Z(l4l)  - 2(1))  - GAMMA ( I ) 4 T 

T • -(TI  - SORT ( TI  • TI  4 T • TEMP))  / T 

CB( I ) • 54575.05416  • T / BLPK(I)  / (1.  4 T • T) 

CBI ( I ) ■ CB( I ) / T 
GO  TO  37 

••SPECIAL  CASE,  GRADIENT  REAL  NUMBER 

TEMP  ■ C(I)  - 2.  • GAMMA ( I ) • (2(141)  - 2(1)) 

TEM  ■ TEMP**2  4 CI(I)**2 
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171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 
227 


XRE  ■ (T  * TEMP  + T I * CI(J))  / TEM 
XIM  - -(T  • CI(I)  “ TI  * TEMP)  / TEM 
TEM  > XRE**2  ♦ XIM**2 
CB( I ) * SORT ( ( XRE  ♦ SQRT(TEM))  • .5) 

CBI ( I ) * .5  * XIM  / CB( I ) 

GAMMAI( I)  « 0. 

GO  TO  33 

29  TEMP  * C(I)  - 2.  * GAMMA(I)  • (2(1+1)  - 2(1)) 

CB(I)  = SORT  (T  / TEMP) 

GAMMA!  ( I ) « .5  * (CI(I)  - TI  / CB(I)**2)  / (20+1)  - 2(1)) 

GO  TO  33 
27  N ■ N - 1 
33  CONTINUE 

C COMPUTE  USEFUll  QUANTITIES 
PRINT  58 

58  FORMAT  (7X.6H  RE  M .8X.6H  IM  M t9X,5H  L/KM.8X.6H  RE  C ,8X, 

* 6H  IM  C .5X.12H  RE  C BOTTOM. 4X, 12H  IM  C BOTTOM. 1 0X ,9H  GRADIENT  ) 
OMEGA  * 6.28318530700  * FREQ 

00  40  I * 1 ,N 

TEMP  * C( I )**2  + C I ( I ) **2 

CAY ( ! ) - OMEGA  * C(I)  / TEMP 

CAY  1(1)  » -OMEGA  * Ct ( I ) / TEMP 

CAY  SQO)  * CA  Y ( I ) * *2  - CAYl  ( I ) +*2 

CAY  S0I(1)  * 2.00  * CAY ( I ) * CAY  1(1) 

TEMDR  s -2.  * ( GAMMA ( I ) * CAY  SQ(I)  - GAMMA  1(1)  * CAY  SQI(I)) 
TEMOI  * -2.  « ( GAMMA ( I ) * CAY  SQI(I)  ♦ GAMMA  1(1)  * CAV  SQ(I)) 

G CU(I)  » (TEMDR  « C(I)  + TEMOI  * Cl ( I ) ) / TEMP 

G CUI(I)  » (TEMOI  * C ( I ) - TEMDR  • CIO))  / TEMP 

TEM1  « OCBRT ( -DSQRT ( GAMMA(I)**2  + GAMMAI ( I ) **2)  • 2.*0MEGA»*2) 

TEM1I  * DATAN  ( ABS( GAMMAI ( 1 ) / GAMMA ( I ) ) )/  3. 

CRTG  * TEM1  * DCOS(TEMII) 

CRTGI  > TEM1  • DSIN(TEMII) 

IF  (GAMMA(I)  .LT.  0.)  CRTG  * -CRTG 
IF  (GAMMAI ( I ). LT.  0.)  CRTGI  > -CRTGI 
G(I)  * (C(l)  * CRTG  + C1(I)  * CRTGI)  / TEMP 
GI(I)  = (C(I)  * CRTGI  - CIO)  * CRTG)  / TEMP 
CON(I)  = G(I)  * C(I)  - GI(I)  * CI(I) 

CON(I)  * 0MEGA»*2  / C0N( 1 ) *»2 
XMI  * — G I ( I ) * (20  +1)  - 20)) 

XM  » -GO)  * (2(1  + 1)  - 2(1)) 

DPK(I)  - -8686.00  * CAY  1(1) 

PRINT  30.  XM.  XMI,  DPK(I),  CO).  Cl  ( I ) . CB(I),  CBI  ( I ) 

* .GAMMAI I).  GAMMAI (I) 

G SOI ( I ) * 2.  * G( I)  • GI( I) 

40  G SO(I)  « G( I ) *»2  - GI ( I ) *»2 
C FINO  MOOES 
NXTRA*0 
IJ  F LAG- 0 
NN  ■ NN  + 1 

IF  (K5  .EQ.  1 ) GO  TO  15 
DO  SO  NN  > 1 ,350 

IS  IF  (IJ  FLAG  .EQ.  1)  GO  TO  53 
52  IF  (NXTRA  .GT.  0)  GO  TO  44 

READ  60,  V, VI .STEP, STEPI, NXTRA 
60  FORMAT  (4D10.4.I10) 

IF  (NXTRA  .GE.O)  GO  TO  62 
V • V + VI  * 1.0-10 
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228 

VI  » STEP  ♦ STEPI  • 1.D-10 

229 

GO  lu  Bb 

230 

62 

IF  (V)  142.301,70 

231 

142 

I F( STEP ) 44,44,143 

232 

C SEARCH  FOR  MODE 

233 

143 

SIZE3  * -1. 

► 

234 

SIZE2*0 

235 

IJ  FLAG* 1 

236 

V«-V 

237 

IF  (NXTRA)  55,55,54 

• 

238 

55 

NXTRA  * 20 

239 

54 

XTRA  * NXTRA 

240 

HOP  * (STEP  - V)  / XTRA 

241 

HOPI *0 . 

242 

lF(STEPl.NE.O. ) HOP I«( STEP I-VI ) /XTRA 

243 

DO  47  IJ  * 1 .NXTRA 

244 

CALL  SETUP 

245 

DET  * VEL 

246 

DET I * VELI 

247 

CALL  DETNT(N, VEL, VELI) 

248 

DELTA  * VEL 

249 

OF  | T T - WFI  t 

250 

SIZE  * DELTA*DELTA  ♦ DELT I*DELTI 

251 

PRINT  56.  V.  VI,  SIZE,  VEL,  VELI 

252 

56 

FORMAT  (2F12.3,  3D17.5) 

253 

IF  ( (SIZE2.LT. SIZE3) .AND. (SIZE. GT.SIZE2))  GO  TO  45 

254 

46 

SIZE3*SIZE2 

255 

SIZE2*SIZE 

256 

V * V ♦ HOP 

257 

V I * VI+HOP I 

258 

GO  TO  47 

259 

45 

V » V - HOP 

260 

TEMP  * HOP  / (SIZE  - SIZE2) 

* 

261 

DELT I * TEMP  * (DET  * VELI  - OETI  * VEL) 

262 

TEMP  » .500  • ( SIZE3  - SIZE)  / (SIZE3  ♦ SIZE  - SIZE2  - SIZE2) 

263 

DELTA  * HOP  * TEMP 

264 

I F ( HOP I . EQ . 0 ) GO  TO  49 

265 

V I * V I-HOP I 

266 

DELT  AI * HOPI *TEMP 

267 

GO  TO  49 

268 

47 

CONTINUE 

269 

IJ  FLAG*0 

270 

NXTRA*0 

271 

GO  TO  52 

272 

53 

S'7'  .*-1. 

273 

S,  * ■ 

274 

Gl  . 46 

275 

44 

NXTf.»  * NXTRA  - 1 

276 

V * 3 • (PHASE  V(NN-1)  - PHASE  V(NN-2))  ♦ PHASE  V(NN-3) 

277 

VI  » 3 . * (PHASI  V(NN-1)  - PHASI  V(NN-2))  ♦ PHASI  V(NN-3) 

278 

STEP  * (PHASE  V(NN-1)  - PHASE  V(NN-2))  • .0001 

279 

70 

CALL  SETUP 

280 

CALL  OETNT (N ,DET , DET I ) 

281 

80 

FORMAT  (/,  2D20.11,  4013. 4) 

282 

VEL  * DET 

283 

VELI  - DETI 

284 

DELTA  « STEP 

I 

I 


L 
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342 

COL (3)  • TEM1 

343 

TEMP  > COLO) 

344 

COL ( 4 ) . (TEM1  - TEMP)  • 1.010 

345 

COL (5)  • -NN 

346 

PUNCH  64,  (COL(t),  I ■ 1,5) 

347 

64 

FORMAT  (5110) 

348 

63 

AF( 1 .NN)  « A( 1 ,3) 

349 

AG( 1 .NN)  >0(1 ,3) 

350 

BF ( 1 , NN ) >-A(1,4) 

351 

BG(I.NN)  >-0(1,4) 

352 

PHASE  V t NN ) > V 

353 

PHASI  V ( NN ) > VI 

354 

IF  (K6  .E0.  1 ) GO  TO  73 

355 

74 

PRINT  81,  V.  VI,  DET,  DET I , SIZE,  CNTR 

356 

73 

LL  • N - 1 

357 

IF  (LL-1)  95,96,97 

358 

96 

I • 0 

359 

GO  TO  98 

360 

97 

00  110  d > 2 , LL 

361 

I > J ♦ d - 2 

?R2 

TEVH"  - A(I , 2 ( ♦ AF(  d-1  , NN ) - Q(I,2)*AG(  d-1  , NN) 

363 

TEMNI  > Q( I . 2 ) * AF ( J-1  ,NN)  ♦ A(I.2)*AG(  d-1  , NN) 

364 

TEMOR  • A ( I . 3 ) >A(  1 + 1 .4)  - 0(I.3)*0<  1 + 1 ,4)  - 

365 

1 *(I,4).A(  1+1  .3)  ♦ 0(1. 4). 0(  1*1  .3) 

366 

TEMOI  * 0( I . 3 ) • A(  M ,4)  + A( I , 3 ) >Q(  1 + 1 ,4)  - 

367 

1 0(1, 4). A(  1 + 1 .3)  - A(  1,4)*0(  1+1  ,3) 

368 

TEMOEN  * TEMOR* TEMOR  ♦ TEMOI *TEM0I 

369 

TEMRNU  * TEMNR*  TEMOR  ♦ TEMNl *TEM0I 

370 

TEMINU  * TEMNI* TEMOR  - TEMNR»TEMOI 

371 

TEMP  * TEMRNU  / TEMOEN 

372 

TEMPI  > TEMINU  / TEMOEN 

373 

BF(d.NN)  * - ( TEMP*  A ( 1*1  ,4)  - TEMPI*Q(  U1  ,4)) 

374 

BG(d.NN)  » -( TEMPI »A(  1 + ) ,4)  ♦ TEMP>0(  U)  ,4)) 

375 

AG(d.NN)  > TEMPI*A(  Ul  ,3)  ♦ TEMP*0(  1 + 1 ,3) 

376 

110  AF(d.NN)  > TEMP.Al  1+1  ,3)  - TEMPI*0(  1+1  ,3) 

377 

98 

TEMNR  « — ( A ( 1+2,2)  • AF(LL.NN)  -0(11-2.3)  * AG(LL.NN)) 

378 

TEMNI  . - ( 0(  I +2  , 2 ) * AF ( LL.NN)  ♦ A(  1+2  , 2 ) • AG( LL , NN ) ) 

379 

TEMOEN  • 4(  1*2  ,3)*A(  1*2  ,3)  *•  0(  1+2  ,3)*Q(  1*2  ,3) 

380 

TEMRNU  • TEMNR* A(  1 + 2 ,3)  ♦ TEMNI>0(  1+2  ,3) 

381 

TEMINU  • TEMN I • A ( 1*2  ,3)  - TEMNR*0(  1+2  ,3) 

382 

BF(N.NN)  > TEMRNU  / TEMOEN 

383 

BG(N.NN)  > TEMINU  / TEMOEN 

384 

95 

AF(N.NN)  > 0. 

385 

AG ( N , NN ) > 0. 

386 

C 

FIND  NORMALIZING  FACTOR 

387 

0( NN ) « 2.1242929600  * RHO(1)>>3  / G( 1 ) 

388 

DI (NN)  * 0. 

389 

DO  111  I > 2 , N 

390 

TEMRSP  * AF ( 1-1  , NN ) >B(  2*1-2  ,2)  - AG(  1-1  ,NN)*BI( 

2*1-2 

.2) 

391 

1 BF ( 1-1  , NN ) • B(  2*1-2  ,1)  - BG(  1-1  ,NN)*Bt(  2*1-2 

.1} 

392 

TEMISP  * AG(  I-t  , NN ) «B(  2*1-2  .2)  + AF(  1-1  ,NN)*BI( 

2*1-2 

.2) 

393 

1 BG(  1-1  , NN) *B(  2*1-2  ,1)  ♦ BF ( 1-1  ,NN)*BI(  2*1-2 

.1) 

394 

AX1  « TEMRSP*TEMRSP  - TEMISP*TEMISP 

395 

AX1 I * TEMRSP  * TEMISP 

396 

AX1 I • AX1I  ♦ AX1 I 

397 

TEMOR  ■ (G(I-1)«*2  + G I ( 1—1 )**2) 

398 

TEMOI  * G ( I ) * * 2 + Gt ( I )**2 

399 

TEMP  - 

( RHO( 1-1  ) 

/ 

RHO(D)  / TEMDI 

400 

TtMl 

■(Z8  (1-1) 

• 

G ( 1—1 ) ♦ ZBl ( 1- 

1) 

• Gl(l-IJ) 

/ TEMDR 

401 

* -(ZT 

(I)  • G(I) 

♦ 

ZTI(J)  • Cl ( I ) ) 

• 

TEMP 

402 

TEM1  I 

■(ZBI(I-I) 

* 

G ( I — 1 ) - ZB  (I- 

D 

• GI(I-I)) 

/ TEMDR 

403 

• -(ZTI(I)  » G( I ) 

- 

ZT  (I)  • 01(D) 

• 

TEMP 

404 

TEMRSP 

• AF(  1-1 

• 

NN)*B(  2*1-1  ,2) 

- 

AG(  1-1  , NN ) * B I ( 

2*1-1 

.2) 

405 

1 BF  ( 

1-1  ■ NN ) * 

B< 

2*1-1  .1)  - BG( 

1 

-1  , NN ) * B I ( 

2*1-1 

.1) 

406 

TEMISP 

• AG(  1-1 

• 

NN) *B(  2*1-1  ,2) 

♦ 

AF(  1-1  , NN) *BI ( 

2*1-1 

.2) 

407 

1 BG( 

1-1  ,NN)» 

B( 

2*1-1  ,1)  ♦ BF( 

I* 

-1  , NN ) * B I ( 

2*1-1 

.1) 

408 

409 

410 

411 
419 

413 

414 

415 

416 
41? 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 

431 

432 

433 

434 

435 

436 

437 

438 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453 

454 

455 


111 


TEMRSP*  TEMRSP  - TEMISP*TEMISP 
TEMRSP  * TEMISP 
AX2I  + AX2I 

/ (G  CU( 1-1 )**2  ♦ G CUI (1-1 )**2> 
RH0( I ) / (G  CU(I  )**2  ♦ G CUI ( I )**2) 

G CU(I-I)  • TEMDR  - G CU(I)  • TEM0I 
G CUI (I)  * TEMDI  - G CUI (1-1)  • TEMDR 
AX1 *TEM1  - AX  1 I • T EMI  I 
AX  1 1 • TEM1  + AX1 *TEM1 I 
AX2  • TEM2  - AX2I  » TEM2I 
AX2I  * TEM2  ♦ AX2  * TEM2I 
D(NN)  ♦ TEMR1  / RHO(I-I)  ♦ TEMR2 
"T'NN)  ♦ TEMI1  / RHO(I-I)  ♦ TEMI2 


AX2 
AX2I 
AX2I 

TEMDR  « RHO( 1-1 ) 
TEMDI 
TEM2 
TEM2I 
TEMR1 
TEMI1 
TEMR2 
TEMI2 
D(NN) 
ni(NNi  - 


112 


130 


CONTINUE 

IF  ( K 1 .GT.  3)  DA(NN)  ■ DSQRT ( ( D(NN) **2  ♦ DI(NN)**2)  * FREQ  / 
• PHASE  *(NN)) 

EIGEN(NN)  « LAMBDA 
El GEN I ( NN)  ■ LAMBDI 
IF  (K6  .GT.  2)  GO  TO  131 
L * 0 
K ■ 24 

DO  112  I » 1 ,N 
L « L ♦ 1 

COL(L)  • SNGL ( 2T( I ) ) • 100. 

L « L ♦ 1 

COL( l)  • SNGL ( 2T 1(1))  * 1000. 

K « K ♦ 1 

COL(K)  • SNGL ( 2B( I ) ) * 100. 

K » K ♦ 1 

COL(K)  » SNGL(ZBKD)  * 1000. 

CONTINUE 


132 

131 


109 


50 


PRINT  130, 
PRINT  130, 
FORMAT  (4H 
M * N ♦ N 
PRINT  132, 


(COL(I),  1-1 . L ) 

( COL ( I ) , 1-25. K) 
2 - . 11(18,15)) 


( RAT  10(1),  I - 1.M) 

FORMAT  ( 1 1 ( 1 X , 2F5. 3 ) ) 

CB  LOSS(NN)  » - LAMBDI  * 8686.00 

PHINV  * V • PHASE  V ( NN-1 ) /((V  - PHASE  V(NN-I))*  FREO) 

PRINT  109.  NN.EIGEN(NN) , EIGEN 1 (NN),D(NN),DI(NN) .PHINV, OB  LOSS(NN) 
FORMAT  ( 3H  N-.I5.10H  LAMBDA  >,2E15.7,4H  D-  .2E15.7, 

• 12H  INT  RANGE  « , F8.0,  6H  L/K  F8.S) 

IF  (K3  . EQ.  0)  GO  TO  50 
CALL  RCOEF  (K3) 

CONTINUE 

C READ  IN  SOURCE  AND  RECEIVEPS  DEPTHS 
301  NRT  ■ NR 
NR  - 0 


51 


456 

NN  • NN  - 1 

457 

MP1  • HI  ♦ 1 

458 

IF  ( K 1 .NE.  3)  GO  TO  321 

459 

NR  - NR  T 

460 

GO  TO  501 

461 

321 

REAO  20,  SOURCE 

462 

320 

NR  • NR  + 1 

463 

REAO  20,  RECVRS(NR),  FINAL,  STEPP 

464 

IF  (NR.GT.50)  GO  TO  300 

465 

350 

IF  (RECVRS(NR)  . EO . 0 . ) GO  TO  300 

466 

310 

IF  (FINAL  .EQ.O. ) GO  TO  320 

467 

330 

RECVRS( NR+ 1 ) « RECVRS(NR)  ♦ STEPP 

468 

IF  (RECVRS(NR-M)  . GT.  FINAL)  GO  TO  320 

469 

340 

NR  • NR  + 1 

470 

I F ( NR  .GT.  50)  GO  TO  300 

471 

GO  TO  330 

472 

300 

PRINT  303 

473 

303 

FORMAT  ( /21H  SOURCE  AND  RECEIVERS  ) 

474 

PRINT  21 .(DEPTH(I) ,1-1 , NR ) 

475 

21 

FORMAT  (8F10.2) 

476 

C COMPUTE  DEPTH  FUNCTIONS 

477 

00  500  I - 1 , NN 

478 

LOC  * 1 

479 

DO  305  J • 1 , NR 

460 

IF  ((J  .EO.  1)  .AND.  (HI  .GT.  5))  GO  TO  305 

481 

LCTR  « 0 

482 

380 

IF((DEPTH(d)  .GE.  Z ( LOC ) ) . AND. ( DEPTH( J)  .LT.  J 

483 

371 

IF  (LOC  .GE.  N)  GO  TO  385 

484 

370 

LOC  * LOC  ♦ 1 

485 

GO  TO  380 

486 

385 

IF  (DEPTH(J)  .GE.  Z ( LOC ) ) GO  TO  360 

487 

390 

LOC*  1 

488 

lctr*lctr+i 

489 

IF  (LCTR  .GT.  2)  GO  TO  305 

490 

GO  TO  380 

491 

360 

XI  - CAY  (LOC)  - EIGEN  (I ) 

492 

X2  ■ CAY  (LOC)  ♦ EIGEN  (I) 

493 

X3  * CAYI(LOC)  - EIGENI(I) 

494 

X4  - CAY  I ( LOC ) + E I GEN I ( I ) 

495 

TEMP  * XI  * X2  - X3  • X4 

4«»4 

TEMP!  - XI  * X4  + X3  * X2 

497 

TEMOEN  * G SOI  LOC ) **2  G SOI(LOC)**2 

498 

ZE  * (TEMP  . GSO(LOC)  TEMPI  • G SQl(LOC))  / 

499 

ZEI  « (TEMPI  * GSQ(LuD  «■  TEMP  « GSQI(LOC))  / 

500 

TEM1  * ZE 

SOI 

IF  (ZE  .GT.  -7.5)  GO  TO  436 

502 

S - CAY (LOC) 

503 

T - CAY  I (LOC) 

504 

DO  437  K ■ 1 ,20 

505 

TEMP  - S**2  ♦ T**2 

508 

TEMPI  * (EIGENI(I)  • S - EIGEN(I)  * T)  / TEMP 

507 

TEMP  - (EIGEN(I)  • S + EIGENI ( I ) * T)  / TEMP 

508 

ZE  ■ ((1.  ♦ TEMP)  • (1.  - TEMP)  ♦ TEMPI—2)  • 

509 

ZEI  - -2.  • TEMPI  * TEMP  • CON(LOC) 

510 

ZR  • ZE  / -7.5 

Bit 

IF  (DABS(ZR-1 . ) .LT.  1.0-3)  GO  TO  438 

• 12 

S - EIGEN(I)  ♦ (S  - EIGEN(I))  / ZR 

TEMOEN 

TEMOEN 


52 


r 


513 

437 

CONTINUE 

514 

438 

IF  (G(IOC)  . LT.  0.)  GO  TO  439 

515 

2E  - G(LOC)  * (DEPTH(d)  - Z(LOC))  ♦ TEM1 

516 

IF  ( ZE  .GT.  -7.5)  GO  TO  442 

517 

F(d)  ■ 1 .0-12 

518 

F I ( d ) • 0. 

519 

GO  TO  305 

520 

439 

ZE  ■ G(LOC)  * (DEPTH(d)  - Z(IOC))  ♦ ZE 

521 

IF  (ZE  .GT.  SLIM)  GO  TO  442 

522 

F ( J ) ■ 1 .0-12 

523 

F I ( d ) • 0. 

524 

GO  TO  305 

525 

442 

ZEI  »GI ( LOC ) • (DEPTH(d)  - Z(LOC))  ♦ ZEI 

526 

302 

CALL  HANKEL(ZE,ZE1  , 1 ) 

527 

F ( d ) * ( AF ( LOC , I ) *H  1 R - AG( LOC , I ) *H1 I ♦ BF( LOC , I )*H2R  - 

528 

1 ♦H2I ) • RHO(LOC) 

529 

F I ( d ) * ( AG( LOC , I ) »H1 R ♦ AF(L0C,I)*H1I  ♦ BG( LOC , I ) *M2R 

530 

1 *H2l)  • RHO(LOC) 

531 

305 

CON  INUE 

532 

IF  (K1  .EQ.  2)  GO  TO  451 

GO  TO  432 

534 

451 

PRINT  270,  DEPTH(NR) 

535 

270 

FORMAT ( 7H 1 DE p TH  , F 5 . 1 , 6X , 3ME-8 , 1 7X , 3HE-6 , 1 7X , 3HE-4 , 1 7* 

536 

* 17X.3HE  0 ) 

537 

432 

IF  (HI  . LT.  4)  GO  TO  431 

538 

IF  (K1  .GT.  5)  GO  TO  433 

539 

SRES(I)  • ( F ( 1 ) « * 2 ♦ F I ( 1 ) **2)  / DA(  I ) 

540 

GO  TO  500 

541 

431 

TEMOEN  « 0( I ) *D( I ) ♦ 01(1 ) *DI ( I ) 

542 

TEMRE  • F(1)»0(I)  ♦ FI ( 1 ) *01 ( I ) 

543 

FO  * TEMRE/TEMDEN 

544 

FDI  » (0(1)  * F I ( 1 ) - DI(I)  * F( 1 ) ) / TEMOEN 

545 

433 

00  400  K • 2, NR 

546 

d « K - 1 

547 

L • d • NN  - NN  + I 

548 

IF  ( K 1 .LT.  6)  GO  TO  435 

549 

FF  . SRES(I)  ♦ (F(K)**2  ♦ FI(K)**2)  / DA(I) 

550 

GO  TO  436 

551 

435 

FF  - FO  • F (K)  - FOI  « FI(K) 

552 

FFI  « FD  • F I ( K ) ♦ FDI  • F (K ) 

553 

436 

UU( L ) * FF 

554 

UUI(L)  * FFI 

555 

452 

GO  TO  (400,410,420,400,400,400,400,400),  X1P1 

556 

C PLOT  DEPTH  FUNCTIONS 

557 

420 

DO  210  II  « 1,120 

558 

210 

COL ( II)-  1 H 

559 

DO  220  11*  20,100,20 

560 

220 

COL(II)«  1 H I 

561 

FE  » FF  * FF  ♦ FFI  * FFI 

562 

IF  ((FE.GT. IE-20). AND. (FE.LT. 10000.))  GO  TO  240 

563 

GO  TO  250 

564 

240 

INT  ■ 100.00  ♦ 2.1714700  * OLOG(FE) 

565 

COL(INT)  ■ 1 H* 

566 

GO  TO  225 

567 

250 

C0L(2)  * 1 H* 

568 

225 

PRINT  260,  COL 

569 

260 

FORMAT  ( 120A1 ) 

53 


570 

071 

572 

573 

574 

575 

576 

577 

578 

579 

580 

581 

582 

583 

584 

585 

586 

587 

588 

589 

590 

591 

592 

593 

594 

595 

596 
5S7 

598 

599 

600 
601 
602 

603 

604 

605 

606 

607 

608 

609 

610 
611 
612 

613 

614 

615 

616 

617 

618 

619 

620 
621 
622 

623 

624 

625 

626 


?R 

GO  TO  400 

INT  DcPTH  FUNCTIONS 

410 

F MAG(J)  « SQRT  (FF  • 

FF  ♦ FFI  • FFI) 

430 

IF  (FF)  430,440,450 

F ANG(J)  - AT AN( FF I / 

FF)  * 57.2957795100  ♦ 

440 

GO  TO  400 

F ANG(J)  « 90. 

450 

GO  TO  400 

F ANG(J)  . AT  AN ( FF I / 

FF)  * 57 . 29577951  DO 

170 

FORMAT  ( 1 0F  12. 4) 

180 

FORMAT ( / 1 0E1 2 . 3) 

400 

CONTINUE 

441 

IF  (K1.EQ.1)  GO  TO  441 
GO  TO  500 

PRINT  180,  (F  MAG(K), 

K . 1,d) 

PRINT  170,  (F  ANG(K), 

K • 1,0) 

500 

CONTINUE 

C CALCULATE  ATTENUATION  AND  READ  IN  RANGES 

IF  ( ( K1  .EG.  4 ) . OR . ( K 1 .EQ.  5))  PRINT  180,  (SRES(K), 
IF  (HI  .EQ.  5)  PUNCH  434,  (SRES(K),  K ■ 1,NN) 

NR.NR-1 

IF  IK?  IT.  31  GO  TO  501 
IF  (K2  .EQ.  4)  K8  - 3 
IF  (K2  .EQ.  3)  KB  « 2 
K2  « 0 

501  KX  » K2  ♦ 1 

GO  TO  (561 .551 .551 ) , KX 

551  PRINT  533.  NN.N,  C(1),  Z(2),  C<2).  Z(3),  C(3>,  2(4), 
* SOURCE,  RECVRS(40),  FREQ 
533  FORMAT  (1H1,  215,  10F10.4) 

ICTR=0 

R LOS ( 1 ) ■ 120. 

LEVEL! 1)  • 1 
DO  562  I « 1,5 
P LEV) I )=40. 

J COU ( I ) =4 
J COUNT ( I ) *-6 

562  CONTINUE 

IF((K2  .EQ.  2). AND. (NR  .GT.  5))G0  TO  772 
GO  TO  561 
772  NR  * 5 
561  NL  « NN 

PRINT  524,  NL 

524  FORMAT  (18.  1 3H  MODES  IN  SUM  ) 

LL  * 1 

IF  ( K9  .GT.  0)  NL  ■ K9 

READ  20,  RANGE,  FINAL  R,  STEP  R 

IF  ( K8  .EQ.  3)  PUNCH  30,  RANGE,  FINAL  R,  STEP  R 

IF  (RANGE)  563,1,564 

563  NN  * NN  + 1 

READ  11,  K1,  K2 , K3 , K4,  K5,  K6,  K7,  KB,  K9 
PRINT  11,  K1 , K2 , K3 , K4,  K5,  K6,  K7,  K8,  K9 
GO  TO  301 

564  IF  (FINALR  . LE.  0.)  GO  TO  550 
FINAL  R . FINAL  R ♦ 1 .0-3 

560  IF  (RANGE  .GE.  FINALR)  GO  TO  561 
550  R ATTEN  « RANGE  • ATTEN  - 9.9429946 


K«  1 , NN) 


C(4), 


54 


627 

IF  ( K 1 .GT.  5)  RATTEN  ■ 0.00 

628 

IF  (K i . ui  . 5)  GO  TO  536 

629 

IF  (K7  .EO.  2)  RATTEN  « 4.342944800  * DLOG(FREQ) 

630 

IF  (RANGE  • DB  LOSS(NL)  .IT.  15.04)  GO  TO  522 

63t 

521 

NL  - NL  - 1 

632 

522 

DO  520  I * 1 , NL 

633 

IF  ( K7  .IT.  2)  GO  TO  523 

634 

FMAG(I)  ■ PHASE  V( I) 

635 

GO  TO  520 

636 

523 

TEMP  RE  * EIGEN(I)  • RANGE 

637 

TEMPIM  * EIGENI ( I ) ‘RANGE 

636 

CALL  HZERO(TEMPRE, TEMPIM) 

639 

HZER02( I ) » HZEROR 

640 

IF  (K7  .EO.  0)  GO  TO  520 

641 

FMAG(I)  - HZER0R**2  4 HZEROI **2 

642 

520 

HZER2I  ( I ) ■ HZEROI 

643 

536 

L « 0 

644 

DO  540  J « 1 , NR 

645 

FF  « 0 

646 

FFI  - 0 

647 

TEMP  * 0 . DO 

648 

DO  530  I * 1 .NL 

649 

K «'L  * I 

650 

IF  ( K1  .LT.  6)  GO  TO  537 

651 

TEMP  * TEMP  4 UU(K) 

652 

GO  TO  530 

653 

537 

IF  (K7  .EO.  0)  GO  TO  534 

654 

TEMP  . TEMP  4 FMAG(I)  • (UUI(K)**2  ♦ UU(K)**2) 

655 

GO  TO  530 

656 

534 

TEMIM  = UUI(K)  * HZER02(I)  4 UU(K)  • HZER2I ( I) 

657 

TEMRE  * UU ( K ) * HZER02( I ) - UUI(K)  * HZER2 1(1) 

658 

FF  * FF  4 TEMRE 

659 

FFI  « FFI  4 TEMIM 

660 

530 

CONTINUE 

661 

IF  (Kt  .GT.  5)  GO  TO  535 

662 

IF  ( K7  .GT.  0)  GO  TO  535 

663 

TEMP  « F F • *2  4 FFI **2 

664 

535 

T RE  * TEMP 

665 

RX  « -4.3429448  * A LOG  ( T RE)  4 R ATTEN 

666 

R LOSS(J)  » RX 

667 

»F  fK4  ,lT.  2)  GO  TO  545 

668 

T RE  > -4.3429448  • ALOG( UU( K ) » *2  4 UUI(K)*»2) 

669 

PRINT  170,  RECVRS(d),  RX.  T RE 

670 

IF  ( K4  .NE.  3)  GO  TO  545 

671 

545 

CONTINUE 

672 

L • L 4 NN 

673 

IF  ( K8  .LT.  2)  GO  TO  540 

674 

IF  (K8  .EO.  3)  GO  TO  538 

675 

LPCH  * -RLOSS(J)  * 10.00  4 1400.500 

676 

IF  (LPCH  .LT.  0)  LPCH  - 0 

677 

IF  (LPCH  .GT.  999)  LPCH  « 999 

678 

LOSPCH(d.LL)  - LPCH 

879 

I RNG  > RANGE  / 1000.00 

680 

IF  (LL  .EO.  25)  PUNCH  903,  IRNG,  ( LOSPCH(d , LLL) , 

681 

903 

FORMAT  (15,2513) 

682 

GO  TO  540 

683 

538 

CONTINUE 

1.25) 


55 


684 

LOSS(d)  • (140.05  - RX)  * 10. 

685 

IF  (LOSS(d)  . LT.  0)  LOSS(d)  • 0 

686 

IF  ( LOSS ( J ) .GT.  989)  LOSS(J)  ■ 999 

687 

540 

CONTINUE 

688 

GO  TO  (770, 780. 716). KX 

689 

C PLOT  DB  LOSSES 

690 

712 

COL( 15) * 1HI 

691 

COL ( 39 ) » 1HI 

692 

COL ( 63 ) * 1HI 

693 

COL ( 87 ) * 1 HI 

694 

COL( 111 ) » 1 HI 

695 

COL ( 27 ) a 1HX 

696 

COL( 51 ) * 1 HX 

697 

COL (75 ) * 1 HX 

698 

COL ( 99 ) « 1HX 

699 

I PLACE  « 135 

700 

DO  787  J > 1 , NR 

701 

I PLACE  « I PLACE  - 24 

702 

783 

I PLOT  . P LEV ( J ) - R LOSS(J) 

703 

IF  (I  PLOT  .GT.  10)  GO  TO  776 

704 

(?*?  TO  777 

705 

776 

P LEV( J | • P LEV(J)  - 20. 

706 

d COUNT ( J ) « d COUNT ( d ) - 2 

707 

d COU(d)«d  C0U(d)-2 

708 

786 

COL(  I PLACE  ♦ 1)  « 1H0 

709 

IF  (P  LEV(d)  - 100.)  778,779,781 

710 

778 

dC  • d COU(d)  ♦ 1 

711 

COL ( I PL ACE ) • ING(dC) 

712 

GO  TO  703 

713 

779 

COL( I PLACE)  * 1 HO 

714 

GO  TO  782 

715 

781 

dC  * d COUNT (d)  ♦ 1 

716 

COL ( I PL ACE ) ■ ING(dC) 

717 

782 

COLO  PLACE  - 1)  « 1 HI 

718 

GO  TO  783 

719 

777 

IF  (I  PLOT  .LT.  -9)  GO  TO  784 

720 

GO  TO  785 

721 

784 

P LEV(d)  « P LEV(d)  + 20. 

722 

d COU(d)*d  COU( d )+2 

723 

d COUNT (d)  « d COUNT (d)  * 2 

724 

GO  TO  786 

725 

785 

IPP  » I PLACE  ♦ I PLOT 

726 

COL(IPP)  * 1H+ 

727 

787 

CONTINUE 

728 

GO  TO  750 

729 

C CONTOUR  LOSS  FIELD 

730 

780 

DO  590  dd  • 1,120 

731 

590 

COL(dd)  * 1H 

732 

COL ( 31 ) « 1HI 

733 

COL ( 61 ) * 1HI 

734 

COL( 91 ) * 1HI 

735 

DO  640  dd  ■ 2,41 

736 

LEV  > 1 

737 

620 

IF  ( RLOS  (dd)  .LT.  CONTR( LEV) ) GO  TO 

738 

GO  TO  610 

739 

600 

LEV-LEV+1 

740 

GO  TO  620 

I 


741 

742 

743 

744 

745 
748 

747 

748 

749 

750 

751 

752 

753 

754 

755 

756 

757 

758 

759 

760 

761 

762 

763 

764 

765 

766 

767 

768 

769 

770 

771 

772 

773 

774 

775 

776 

777 

778 

779 

780 

781 

782 

783 

784 

785 

786 

787 

788 

789 

790 

791 

792 

793 

794 

795 

796 

797 


610  IF  (LEV  .EQ.  LEVEL(JJ))  GO  TO  640 
650  IF  (LEV  .01 . LEVEL ( JJ) ) GO  TO  660 
GO  TO  670 

660  II  • LEVEL(JJ) 

GO  TO  680 
670  11  « LEV 

680  JJJ  « 124  - 3* JJ 

COL(JJJ)  • J SMBL ( II) 

LEVEL! JJ > » LEV 
640  CONTINUE 

COL( 1 ) « 1 HI 

PRINT  261,  (COL(II),  II  ■ 1.119) 

716  00  690  JJ  > 1,120 

690  COL(.tJ)  ■ 1H 

ICTR  x ICTR  ♦ 1 

IF  (ICTR  .EQ.  10)  GO  TO  700 

GO  TO  714 

700  TEMP  x (RANGE  ♦ 1 . ) / 10000. 

INO  » TEMP 

COL( 2 ) < ING( IND+1 ) 

TEMPI  • IND 

TFMD  , (TfMO  - TFMP1 ) * 10. 

INO  < TEMP 

COLO)  x 1NG(  IND+1) 

TEMPI  . IND 

IND  > ( TEMP  - TEMPI ) • 10. 

COLO)  x ING(  IND+1) 

C0L(4)»1H. 

C0L(6)x 1HK 
COL ( 7 ) « 1 HV 
C0L(8)> 1H0 
COLO)  • 1HS 
ICTR«0 

714  GO  TO  (710, 712), K2 
710  COLOt  ) * 1HI 
COL ( 6 1 ) * 1 HI 
COL ( 91 ) * 1HI 
DO  720  JJ  < 1,40 
TEMP  « LEVEL ( JJ) 

TEMPI  x o. 

830  'F  ( LEVEL ( JJ ) .GT.  LEVEL(JJ+1))  GO  TO  730 
> TO  740 

730  II  x LEVEL! JJ ) - 1 
KK  x 1 

860  EX  > (CONTR(II)  - R LOS  (JJ)  )/  (RLOS  (JJ+1  ) - CONTR(II)) 
DO  760  LL  « 1,3 

IF  (EX  .LT.  TEST ( L L ) ) GO  TO  800 
760  CONTINUE 
LL  < 4 

800  JJLL  • 125  - 3* JJ  - LL 
COL(JJLL)  > J SMBL (It) 

GO  TO  (810, 820), KK 
810  LEVEL(JJ)  - LEVEL(JJ)  * 1 
GO  TO  830 

740  IF  (LEVEL(JJ)  .LT.  LEVEL(JJ+1))  GO  TO  840 
GO  TO  720 

840  II  ■ LEVEL! JJ) 
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1 


798 

KK«2 

799 

uu  iu  bt>u 

800 

820 

LEVEL(Jd)  • LEVEL! dd)  ♦ 1 

801 

GO  TO  740 

802 

720 

LEVEL(dd)  * TEMP 

803 

COL( 1 ) » 1 HI 

804 

750 

PRINT  261,  (COL(ll),  11  • 1.119) 

805 

261 

FORMAT  (IX,  119A1) 

806 

GO  TO  581 

807 

C PRINT  DB  LOSSES 

808 

770 

PRINT  580,  RANGE,  (R  LOSS(K),  K • 1 ,NR) 

809 

LL  ■ LL  ♦ 1 

810 

IF  (LL  .GT.  25)  LL  - 1 

811 

580 

FORMAT  (F9.0,  2X,  18F6.1) 

812 

581 

RANGE  * RANGE  ♦ STEP  R 

813 

IF  (K8  ,NE.  3)  GO  TO  560 

814 

PUNCH  980,  (LOSS(I),  l»  1.NR) 

815 

980 

FORMAT  (2613) 

816 

GO  TO  560 

817 

999 

STOP 

818 

END 

58 
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t SUBROUTINE  SETUP 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  LAMBDA . LAMBDI 

COMMON  /HAN/  H2R . H2 1 , HI R , H 1 1 ,H2PR , H2PI ,H1 PR , HI  PI , EXPONT 
COMMON  /EXPO/  EXSUM,  CNTR,  RAT  10(25) 

COMMON/DETMNT/  A(25, 4) ,0(25,4) 

COMMON/ INPUT/  Z( 12 ) ■ N,  OMEGA,  V,  VI,  C0N(12),  GS0(12), 
1 CAY (12),  LAMBDA,  LAMBOI , G(12) 

2 , RHO( 12),  G I ( 12),  G SOI (12),  CAYI(12) 

COMMON  /LIMIT/  TLlM,  EXPON  , SLIM 
COMMON/PARTS/  ZT ( 1 2 ) . ZTI ( 1 2 ) , ZB( 1 2) , ZBI ( 1 2 ) 

DENOM  . V * V ♦ VI  * VI 
LAMBDA  > OMEGA  * V / DENOM 
LAMBOI  > -OMEGA  • VI/  DENOM 
M ■ N - 1 
DO  10  I « O.M 
IF  (I  .EO.  0)  GO  TO  35 
IF  (ZR  .GT.  -7.4)  GO  TO  25 
IF  (G( I ) ,LT.  0. ) GO  TO  25 
ZE  • G(!)  • (Z(t+1  ) - Z(I ) ) + ZE 
IF  (ZE  .LT.  -7.5)  ZE  ■ -7.5 
GO  TO  26 

25  CONTINUE 

ZE  • G< I)  • (2(1*1  ) - Z(I))  * ZR 
IF  (ZE  .LT.  SLIM)  ZE  - SLIM 

26  CONTINUE 

ZO  • GI(I)  * (Z(f+1)  - Z(  I ) ) * ZI 
30  ZB( I ) < ZE 
ZBI ( 1 ) * ZO 
CALL  HANKEL(ZE.ZQ.O) 

ZB( I ) • ZE 
ZBI ( I ) * ZO 
RAT I0( 2*1)  » EXPONT 
A( 2* I , 1 ) . H2R  * RHO(I) 

0(2*1 , 1 ) * H2I  • RHO( I ) 

A( 2* I , 2 ) * HI R * RHO(l) 

0(2*1,21  * HI  I * RHO( I ) 

A ( 2*  I + 1 , 1 ) > H2PR  * G ( I ) - H2PI  • G I ( I ) 

0(  2* 1 + 1 , 1 ) « H2PI  • G ( I ) + H2PR  • G 1 ( 1 ) 

A ( 2* 1 + 1 , 2 ) » HI  PR  • G( I ) - HI  PI  • GI(I) 

012*1  + 1.2)  - HI  PI  * G( I ) + HI  PR  • GI(I) 

35  CONTINUE 

GSABS  * G SO ( I + 1 ) • • 2 + G SQI(1+1)**2 
XI  « CAY (1+1 ) - LAMBDA 
X2  ■ CAY ( 1+1 ) + LAMBDA 
X3  -CAY  I (1  + 1 ) - LAMBDI 
X4  -CAY  1(1  + 1)  + LAMBDI 
X > XI  * X2  - X3  * X4 
Y • XI  * X4  + X3  • X2 

ZT ( 1 + 1 ) - (X  • G S0(I  + 1)  + Y • G SOU  1+1 ))  / GSABS 
ZT I ( 1 + 1 ) > (Y  * G S0(I+1)  - X * G SOI (1  + 1))  / GSABS 
ZR  - ZT ( 1+1 ) 

ZI  - ZTI ( 1+1 ) 


ZE  - ZR 
56  ZO  • ZI 
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57  IF  (ZR  .GT.  -7.5)  GO  TO  40 

S « CAv ( i + 1 ) 

T ■ CAY  1(1  + 1 ) 

CON  > ( G( I +1 ) • S ♦ GI ( 1+1 ) • T)  / (S**2  ♦ T**2) 

CON  « 1 . / C0N**2 
DO  36  J » 1,20 
TEMP  * S**2  + T**2 

TEMPI  * (LAM8DI  * S - LAMBDA  • T)  / TEMP 
TEMP  • (LAMBDA  * S ♦ LAMBDI  * T)  / TEMP 
ZR  * ((1.  + TEMP)  * ( 1.  - TEMP)  ♦ TEMPI**2)  • C0N(I+1) 

R * ZR  / -7.5 

IF  ( DABS( R-1 . ) .LT.  1.D-3)  GO  TO  41 
S » LAMBDA  + (S  - LAMBOA)  / R 
36  CONTINUE 

41  ZI  * -2.  • TEMPI  * TEMP  * C0N(I+1) 

ZT( 1+1 ) « ZR 
ZTMI  + 1 ) • ZI 


40  CONTINUE 

CALL  HANKEL( ZR , ZI , 0 ) 

ZT ( 1+1 ) * ZR 
7TT ( 1+i > » 71 

RAT  10(2*1  + 1 ) = EXPONT 
IF  (I  .NE.  0)  GO  TO  45 
A( 1 ,3)  = H2R  » RHO ( 1 ) 

0(1,3)  = H2I  * RHO ( 1 ) 

A( 1 ,4)  = HI R * RHO ( 1 ) 

0(1,4)  * HI  I + RHO ( 1 ) 

GO  TO  10 
45  CONTINUE 

A ( 2* I , 3 ) »-H2R  * RHO( 1 + 1 ) 

0(2*1, 3)  *-H2 I * RHO( 1+1 ) 

A( 2* I , 4 ) *-H1R  * RHO( 1+1 ) 

0(2*1. 4)  *-H 1 I * RHO( 1 + 1 ) 

A(2*l+1  ,3)  »-H2PR  * G(  1+1  ) + H2PI  • GI(1*-1) 

0( 2* 1+1 ,3)  *-H2PI  * G( 1+1 ) - H2PR  • GI ( 1+1 ) 

A( 2* 1 + 1 , 4 ) *-H1PR  * G( 1 + 1 ) + HI  PI  • GI ( 1 + 1 ) 

0( 2* 1 + 1 , 4 ) *— HI  PI  * G( 1+1 ) - HI  PR  • GI(I+1) 

10  CONTINUE 

A(  2*N-2  ,4)  « 0. 

Q(  2*N-2  ,4)  ■ 0. 

A(  2*N-1  ,4)  * 0. 

RETURN 

END 


1 


1 

SOMnOu i INt  utTNUN.OET.DETl) 

3 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

9 

DOUBLE  PRECISION  DET,  DET1 

4 

COMMON  /EXPO/  EXSUM,  CNTR,  RATtO(25) 

S 

COMMON/ DETMNT/  A(25,4),  0(35,4) 

• 

OLOSS  « 1. 

r 

CNTR  ■ 0. 

• 

DET  - A( 1 , 3) 

• 

DETI  ■ 0(1.3) 

10 

LIM  • N ♦ N - 3 

11 

DO  100  1-1, LIM, 2 

13 

J • I 

13 

K a 4 

14 

L ■ I 

IS 

M ■ 3 

IS 

II  ■ 1 

17 

GO  TO  500 

IS 

10 

J • I ♦ 1 

19 

K ■ 3 

30 

L • J 

31 

M - 1 

33 

GO  TO  600 

33 

30 

J • I ♦ 2 

34 

L ■ J 

35 

II  ■ 2 

36 

GO  TO  600 

37 

40 

L ■ I ♦ 1 

38 

M ■ 2 

39 

GO  TO  500 

30 

50 

K ■ 3 

31 

M ■ 3 

33 

It  - 3 

33 

GO  TO  600 

34 

60 

K ■ 4 

35 

IF  (I  .EQ.  LIM)  GO  TO  70 

36 

M • 4 

37 

II  ■ 4 

38 

GO  TO  600 

39 

70 

K - 3 

40 

II  * 3 

41 

GO  TO  700 

43 

500 

C - A(L,M)*A(L,M)  ♦ 0(L,M)*Q(L,M) 

43 

80 

8 • <A(J,K)«A(L,M)  ♦ Q(  v) , K ) *Q(  L ■ M) ) / C 

44 

BI  - (O(J.K).A(L.M)  - A(J,K)*0(L,M))/C 

45 

GO  TO  t 10.50),  II 

46 

600 

TO  . A(  <J  , K ) - ( A(  L ,M)*B  - Q(L,M)*BI) 

47 

TDI  ■ 0(J,K)  - ( A( L ,M) *81  ♦ 0(L,M)»B) 

48 

TEM  ■ TO*  *2  TDI  * »2 

49 

TEMP  « A(J,K)**2  ♦ Q( J , K) * *2 

50 

TEMP  » TEM  / TEMP 

51 

IF  (II  .EQ.  2)  GO  TO  92 

S3 

IF  (It  .EQ.  4)  GO  TO  92 

S3 

9(J,K)  ■ Q ( J , K ) * 10.0-18 

54 

A(J,K)  ■ A( J ,K)  « 10.0-18 

55 

IF  (TEMP  .GT.  10.0-35)  GO  TO  92 

56 

CNTR  ■ CNTR  ♦ 1. 

61 


57 

GO  TO  90 

58 

92 

A(d,K)  « TD 

59 

Q(d,K)  » TDI 

60 

90 

GO  TO  (700, 40, 60, 70), II 

61 

700 

C « DET*A( d ,K)  - DETI*Q(d,K) 

62 

OET I • DET *Q( J , K)  ♦ DETI* A( d ,K) 

63 

DET  « C 

64 

GO  TO  (30,100),  II 

65 

100 

CONTINUE 

66 

RETURN 

67 

ENO 

SUBROUTINE  RCOEF  (K3) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON/INPUT/  Z ( 1 2 ) . N,  OMEGA,  V,  VI,  GCUM2),  GS0(12), 
f CAT (12),  LAMBDA,  LAMBDI , G(12) 

2 , RKO( 12),  G I ( 12),  G SO I ( 1 2 ) , CAYI(12) 

DIMENSION  RR( 12) , R I ( 1 2 ) . RA(12),  RT(12),  CYSQ(12),  CYSOI (12) 
COMMON/ RE FL/  AF(12,200),  AG(12,200),  BF(12,200),  BG(12,200), 

2 E I GEN ( 350 ) , E I GEN  1(350), BR (25, 4),  BI(25,4),  CB( 12),  CBI(12), 

3 CAYSQ( 12),  CAYSOI ( 12) , NN 
NM  « N - 1 

I « K3 

IF  (I  . GT.  NM)  I * NM 
110  J » I + I 
K « d ♦ 1 

IF  (NN  .NE.  1 ) GO  TO  102 
L * I ♦ 1 

TEMP  = C B( I ) * • 2 + CBI ( I )• *2 
CY  * OMEGA  • CB( I ) / TEMP 
CYI  * -OMEGA  * CBI (I)  / TEMP 
CVSQ(I)  = CY* *2  - C Y I **2 
CYS'’'1'I'  - CY  * CYI 
CYSOI(I)  * CYSQI(l)  + CYSOI  ( I ) 

102  EL  SO  * C YSO(I)  - E IGEN( NN ) **2  + E IGENI (NN  ) **2 
EL SO  I = C YSQI(I)  - 2. DO  » EIGEN(NN)  * EIGENI (NN) 

TEMP  * ELSQ  + DSORT  (ELS0*»2  ♦ ELS0I**2) 

IF  (TEMP  ,LE.  0 . DO ) GO  TO  107 
EL  * DSORT  (TEMP  * .5D0) 

ELI  * ELSOI  / (EL  ♦ EL) 

103  A * AF( I , NN ) * BR( J , 2 ) - AG ( I , NN ) • BI ( J , 2 ) 

* ♦ BF( I , NN ) • BR ( d , 1 ) - BG( I,NN)*BI(d,1 ) 

B * AF  ( I , NN ) * BI  ( J , 2 ) + AG  ( 1 , NN  ) * BR  ( J , 2 ) 

* ♦ BF(  I , NN ) » B I ( d , 1 ) + BG(  I ,NN)*BR(d,  1 ) 

E * AF ( 1 , NN )*BR(K,2)  - AG( I , NN ) * BI ( K , 2 ) 

* ♦ BF( I , NN ) *BR(K , 1 ) - BG( I ,NN)*BI(K,1  ) 

F * AF(  I , NN ) * B I ( K , 2 ) + AG(  I , NN ) * BR( K , 2 ) 

* «•  BF  ( I , NN  )*BI(K,1  ) + BG(  I ,NN)*BR(K,  1 ) 

C « (F  * EL  - E ♦ ELI)  / (ELSO  * ELSOD 
D « -(E  * EL  ♦ F * ELI)  / (ELSO  + ELSQI) 

TEMP  * (A  + C)**2  ♦ (B  ♦ D) **2 

RR(  I ) » ( A»*2  - C* *2  ♦ B* *2  - D**2)  / TEMP 
R I ( I ) = -2. DO  *(A*D-B*C)/  TEMP 
10  FORMAT  (10D13.5) 

RA(  I ) * 0 

IF  (CB( I ) .GT.  V)  GO  TO  104 

RX  * CB( I ) / V 

RA ( I ) * ACOS(RX)  * 57.296 

104  RT(  I ) * RR( I ) »*2  ♦ R I ( I )**2 
RT ( I ) « 1 .DO  / RT( I) 

R I ( I ) * -DATAN2  ( R I ( I ) . RR ( I ) ) • 57.296D0 
RR( I ) * -4.34294D0  * DLOG  (RT(I)) 

IF  ( K3  . NE . 1 ) GO  TO  108 
I « I ♦ 1 

IF  (I  .LT.  N)  GO  TO  110 

105  CONTINUE 

PRINT  106,  ( RR( I ) , I ■ 1 , NM) 

PRINT  106,  ( RI ( I ) , I ■ 1 , NM) 
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57  PRINT  106.  ( RT( I ) , I ■ 1,NM) 

58  PRINT  106,  ( RA( I ) , I > 1,NM) 

59  106  FORMAT  (9G13.4) 

60  RETURN 

61  108  PRINT  109,  Z(l),  RR(I),  RT(I),  RI(I),  RA(I) 

62  109  FORMAT  (9H  AT  DEPTH,  F7.0.8H  YD,  R ■ ,F9.4  ,7H  DB,  OR, 012. 4, 

63  • 8H,  PH  A ■ , F9.3.16H  DEGREES,  GR  A -.F8.2.9H  DEGREES.  ) 

64  RETURN 

65  107  EL  ■ O.DO 

66  ELI  » DSORT  ( 0A8S  ( ELSQ) ) 

67  GO  TO  103 

68  END 
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SUBROUTINE  HANKEL  (ZR,  Z 1.  IH) 

COMMON  /HAN/  H2R,H2I,H1R,H1l,H2PR,H2PI,H1PR,H1PI,R 

INTEGER  FIPS,  FIQUAD 

DOUBLE  PRECISION  ZMLA2 , TLIM.  R 

COMMON  /LIMIT/  TLIM, EX PONT 

DOUBLE  PRECISION  ZR, ZI , H2R ,H2I ,H1 R , HI  I ,H2PR ,H2PI ,Ht PR ,H1 PI , AO , 

1 A, BO, B, CO, C. DO.O,  C4,  04 , K ,K1 . K2 .C0N4 , STORE  1 . ST0RE2 , 

2 STORE  3, ST0RE4, STORES. ST0RE6 , ST0RE7 . STOREB . ST0RE9 . STOR 1 0 . STOR 1 1 , 

3 ST0RI2,ST0R13,ST0R14,ST0R15.ST0R16,ST0R18,C1 . C2 , C3 ,CPR ,CPI , CTHR 

4 CTP.FR.FI.FPR.FPI  , F t R , F II , F2R , F2 1 ,GR ,G t ,GPR , GPt ,G1 R , 

5 G1 I .G2R.G2I ,H11R,H11l,H12R,H12I,H11PR,H11Pl,H12PR,H12PI, 

6 H21PR.H21PI .H22PR.H22PI .H21R.H21 I .H22R.H22I , PI , SR , SI , SPR . 

7 SPI.S2,  THR,X, XR  , X I ,XPR,XPI ,YR,Yl • ZM, ZMSQ, ZM1 R . ZM1 I , 

8 EXPONT , ABK , ZRTM2R , ZRTM2 I . ZRTM4R . ZRTM4 I , ZRT2R ,ZRT2 I , ZRT2M , ZRT4R 

9 ZRT4I,  Z32R,Z321.STP,ST0R17.STHR,CS,05,FP112 

DIMENSION  A(40) ,B(40)  ,C(40) ,0(40) , C4 ( 20 ) , D4 ( 20 ) ,ZMLA2( 40 ) , 

1  XPR( 40 ) , XP I ( 40 ) , C5( 20) , 05(20),  ZMLA5I20) 

DATA  (A( I), I* 1,36)  / 


t 

2 

3 

4 

5 

6 

7 

8 
9 
A 
B 
C 
0 
E 
F 
G 
H 


■1 . 550 7 2786 15487 157 D-001 . 


-7 . 1 792936553 1 8 1 28300-05 , 
-2.589933497589512370-09, 
-2 . 0 15198799867  345450-  1 4 , 
-5.200459349754700470-20, 
-5.660548752345328790-26, 
-3.031375850066045880-32. 
-8.890812451067134410-39, 
-1 .545475672901 393130-45, 
-1.691724586734780180-52, 
-1.223472353654655730-59, 
-6.078254280234251460-67, 
-2.1423727531 17290340-74, 
-5.504713273139644150-82, 
-1 .055260349669891260-89, 
-1 .53977168218007537D-97, 
-1 .740204228758308300-105, 
-1  .546877903396463700-1  13, 


5. 1 69092871 82905237D-03 , 


5 . 438860344937975960-07 , 

8 . 463834959442850890-1  2 , 
3.650722463527799730-17, 
5.977539482476667200-23. 

4 . 49249900979787999D-29 , 

1 .760380865311292610-35, 

3 . 940962965898552490-42 , 

5 . 399984880857418350-49 . 
4.778883013375085270-56, 
2.851916908285910780-63, 

1 .189016877980096140-70, 
3.567054200994489410-78, 
7.895457936230126430-86, 

1 .317428651273272490-93, 

1 .688346142741310720-101 , 
1 .6891 90670 508938 36D-1  09 , 
1 .338592855137126780-1  1 7/ 


DATA  (B( I ) , 1-1 ,36)  / 


-5 . 65 2489376 2 02298 90- 


• 

1 -1 .495367559841878020-05,  9 

2 -3.994037285902451990-10,  1 

3 -2 . 52 7807 7 0480b 4 93500“ 1 5 , 4 

4 -5.572768308656290780-21,  5 

5 -5.340663090733033160-27,  4. 

6 -2.570196682611954820-33,  1 

7 -6.875088092327653980-40,  2 

8 -1.102217825003022680-46,  3 

9 -1.122556300047279290-53,  3 

A -7.606879255893286000-61,  1. 

B -3.561563187213418130-68,  6 

C -1.188804503195485240-75,  1 

D -2.904623697738803090-83,  4 

E -5. 31 361078500669380D-91  , 6 

F -7 .421557 1 45296762240-99 , 7 

0 -8.050389141952994550-107,  7 

H -6.884688783453903250-1 15,  5 


002,  1 .34583080385769022D-03 

5856894861 65BR4  760-nR . 
167847159620600000-12, 
213012841344155830-18, 
992223987802463200-24, 
009506824874649520-30 . 
4231432351 1 1824370-36, 
923081671908016150-43, 

711 171127956305320-50, 
067093715976172910-57, 
720235019424080960-64, 
7761856682 1 4265850-72 , 
9292510600381 1301D-79, 
068100417001 12477D-87, 
487925256418429550-95, 
969885250533464600-103, 

66265861 598419431D-1 1 1 , 
84835948305632284D-1 19/ 
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M 

DATA  (C( 1) , I • 1 .361  / 

99 

• 

-3. 10145572309743140- 

002,  6.461366089786315470-04 

60 

1 

-6.526633413925571180-06, 

3.884900246384268560-08, 

61 

2 

-1  .523490292699713160-10. 

4.231917479721425440-13, 

62 

3 

-8.76173391 246671 9350- 1 6 . 

1 .404124024433769130-18, 

63 

4 

-1 .793261844743000160-21 , 

1 . 06790 1 08H 2 73958600-24 , 

64 

9 

-1 .6172996435?  ,'236000-27, 

1 . 182236581 525757890-30 , 

69 

6 

-7.393599634307428970-34, 

4.000865602980210480-37, 

66 

7 

-1  .8916622236J1 305190-40, 

7 . B81 92593 17971 04970-44 , 

67 

8 

-2. 9 1 5991  83566.100591 0-47, 

9.642830144368247050-51 , 

60 

9 

-2.867329800026051 160-54, 

7 . 70 7 8 7582 HO 24 331 000-58 , 

69 

A 

-1  .882265159460701120-61 , 

4. 193995453361 6335  ID-65, 

70 

0 

-8 . 5609? 1521 456692200-69 , 

1 . 606 7 79564 8 37967 76D- 7 2 , 

71 

C 

-2. 7823022767  75701 74D- 76 , 

4.45881775124311 1 760-80. 

72 

0 

-6.63218466643 1306210-84, 

9. 180765042128053990-88, 

73 

E 

-1.186685786145945240-91 , 

1 .431 987664427470 1 0D-95 , 

74 

r 

-1  .620812297031658290-99. 

1 .722802186470725220-103, 

73 

Q 

-1 .7229744839191 17130-107, 

1 .624221798566206890-1  1 1 , 

76 

77 

70 

M 

-1 .445680283548096920-115, 

1 .216902595579206160-1  19/ 

data  t • i , ) / 

79 

• 

-2. 260"957504009195D 

-001.  9.420815627003831550-0 

BO 

i 

-1 .495367559841078020-04, 

1 .246139633201565020-06, 

81 

2 

-6. 39045965744  3923180-09, 

2.21 09096032791 39990- 1 1 , 

62 

3 

-5.561 170950574285690-14, 

1 .053253210336038060-16, 

83 

4 

-1.560375126423761420-19, 

1 .867589436218763590-22, 

84 

5 

-1  .815825450849231270-25, 

1 .483517525203620320-28, 

85 

6 

-1 .02807867304  1781930-31 , 

6.  119515910980644810-35, 

86 

7 

-3. 162540522470720830-38, 

1 .432310019234927910-41 , 

87 

9 

-5.731532690015717940-45, 

2.041144 120375967930-48 , 

88 

9 

-6.51 082654027 12 1 9860-52 , 

1 .070927166745465480-55, 

89 

A 

-4.868402723771703040-59. 

1 . 152557463014134240-62, 

90 

B 

-2.493094231049392690-66, 

4.946615537796414070-70, 

91 

C 

-9.034914224285687000-74, 

1 .52410833/430109280-77, 

92 

D 

-2.381791432145810530-81  . 

3 . 467885354450956060-85 , 

93 

E 

-4 .676977490805090540-89 , 

5.904011983340770890-93, 

94 

F 

-6 . 97626371 6578956500-97 , 

7.730788693017460660-101 , 

99 

0 

-8.050389141952994550-105, 

7.892538374463720140-109, 

96 

97 

98 

M 

-7.2977701 10461 137440-113, 

6.374711836531391900-1  17/ 

DATAIC4 ( I ) , 1*1,19)  / 

99 

• 

.1041 6666666666666660000 . 

- . 58763744212962962940000 . 

100 

• 

-.22907160539343377120001 , 

- .51 152469146043830390001  . 

101 

• 

-.900284  76638  740308390001  , 

-.  14134204350396378960002. 

102 

• 

-.203296781761 17332570002, 

-.27649485411 187761090002, 

103 

• 

- . 36093767 1 25929491870002  , 

-.45662621146185479160002, 

104 

• 

-.563561 18497383940990002, 

-.68174312623273330360002 . 

109 

• 

-.811 17244833084638490002, 

-.95184947011821749270002, 

106 

• 

-. 1 10377446 9890957246D003, 

-. 12669480225846890900003, 

107 

• 

-.1441391 3537 1 00938690003 , 

- .16273270737519703760003 . 

108 

• 

-. 1826444261 1 464413830003/ 

109 

110 

DAT  A ( 04 ( 1 ) . 1-1.19)  / 

• -.14583333333333333330000, 

-.52426938657407407340000, 

111 

• 

-.219001074001 06261 22D001 , 

-.49862280807822504170001  , 

<12 

♦ 

-.89102698762513757310001  , 

-.13961075131923849560002, 

<12 

• 

-.201381 05970329288470002 , 

-.27441048801201 1921 1D002, 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 


* -.3586970305443466737D002,  - .45423931 71 1 94671 533D002 . 

* - . 5b1 0J63O448 05 7570000002 , -. 679087439672985 1 5690002 , 

* -. 808391 9802992951 438D002 , -.94894955810820384810002, 

* 1 1007598835585749100003,  - . 1263823305425793323D003 , 

* -.14381582552263232700003,  -.16241253560510208670003, 

* -. 1 82643081 050056686 10003/ 

DATA(C5( I ) , 1*1,19)  / 

* -.80208333333333333220000,  -.22855450236966824530001, 

* -.37786353593898852400001,  -.52746231607110598620001, 

* -. 6771 926 1 70404923857D001 , -. 826995494 1340274075D001 , 

* -.97684336200979706920001,  -.11267213747680068060002, 

* -.12766207438002797210002,  -.14265358922464757740002, 

* -.15764630887783273850002,  -. 17263997303951 79650D002 , 

* -. 18/6343942052581 175D002,  -.20262943441407857240002, 

* -.21 7624 96686 03070326D002,  -.23262048426053788200002, 

* -.24761042850013399850002,  -.26254300552690594930002, 

* -.27720979241646002400002/ 

DAT  A ( DS ( 1 ) , 1*1,19)  / 

* - . 67 708 333333 3 333332 2D000 , -.22029147982062780150001, 

* -.37125903089615039470001,  - .521801 0289498877822D001 . 

* -.67215926158079361730001,  - .8224 1 86533352991 837D001 , 

* -.97261769556783931290001,  - . 1 1 2277b6967 105 /B7810002 , 

« -.12729075204990961410002,  -.14230176248921136010002, 

* -.15731119656707000070002,  -.17231939820272207020002, 

* -.18732661405573092580002,  -.20233302328193519550002, 

* - . 2 173387494^4898291 00002 , -.23234359923945613750002, 

* - . 247340463SB95374440D002 , -. 262225287708803557 1 D002 , 

* -.27586862583465526640002/ 

DATA  (2MLA5(1).  1*1,17)  / 1.E9.715.,  207.,  103.,  47.,  36.4,  27., 

* 22.6,  18.5,  16.6,  14.7,  14.,  12.9,  12.2,  11.5,  10.8,  9.2  / 

DATA  LA2,  LA5  /36, 17/ 

DATA  <2MLA2( 1 ). 1*1 ,40)  / 

1 2. 694430 1 0-1 2, 4. 7348244D-6, 7. 08037 130-4, 9. 63989320-3, 
24.92714940-2,1 .52673010-1 ,3.53247720-1 ,6.78352770-1 , 

31 . 147521500,1 .7731 14100,2.561774900,3.995718100,5.215932700, 
46.602070200,8. 149097200,9.85141 1200,1 .2320405D1 ,1 .491792301 , 

51 . 716363401 , 1 .954030901 ,2.303024601 ,2.644684401 ,2.929282001 , 

63 . 3549744D1 ,3.754124601 ,4.080298001 .4.5784933D1 ,5.028713301 , 

75. 39171  66D1 .5.9582285D1 ,6.453785801 ,7.070137701 ,7.597847201 , 
88.016339901 ,8.694576601 .9.2594255D1 .9.983446D1 , 

91 .057519302,1 .103982802,1 .182016902/ 

DATA  Cl  / 0.57735  02691896257600  / 

DATA  C2  / 0.66666  666666666666D0  / 

DATA  C3  / 0.86602  54037844386400  / 

DATA  PI  / 3.1415926535897932400  / 

DATA  FPI12  / 1.3089969389957471800/ 

DATA  CON4  / .7071067811 865475244D0/ 

AO  « 9.3043671 6929229427D-01 
BO  ■ 6.782987251 44275871D-01 
CO  ■ 4. 652 1 8358464 61 47 14D- 01 
DO  ■ 6.782987251442758710-01 
K - 0 . 85366721 8838951 56D0 

ZMSQ  - ZR*ZR  +ZI*ZI 
RZR  * ZR 
TEMP  - Z1 

IF  (TEMP  .IT.  0.)  TEMP  - -TEMP 
IF  (RZR  .IE.  0.)  GO  TO  51 
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171 

IF  ( RZR  .GT.  4.4)  GO  TO  120 

172 

TEM1  « 7.  - .2632  * RZR**2 

173 

IF  (TEMP  .GT.  TEM1 ) GO  TO  120 

174 

GO  TO  53 

175 

51 

IF  (RZR  .LT.  -9. ) GO  TO  120 

176 

TEM1  » 4.4  + .1375  * RZR 

177 

IF  (TEMP  .GT.  TEM1 ) GO  TO  120 

178 

53 

FIPS  • 1 

179 

ST0RE1  = ZR»ZR-ZI*ZI 

180 

ST0RE2  = 2 . *ZR*ZI 

181 

XR  = STOREWR  -ST0RE2*ZI 

182 

XI  = ST  0~E  1 *Z I +ST0RE2*ZR 

183 

DO  55  MLS=1 . LA2 

184 

IF  (ZMSO  - ZMLA2 (MLS) ) 62,62,55 

185 

55 

CONTINUE 

186 

62 

FR  = AO 

187 

FI  = 0.0 

188 

XPR( 1 ) = XR 

189 

XPI ( 1 ) * XI 

190 

DO  65  M = 1.MLS 

191 

FR=FR+A(M)*XPR(M) 

192 

FI*FI+A(M)*XPI(M) 

193 

XPR ( M+ 1 )=XR»XPR(M)-XI*XPI(M) 

194 

XPI (M+1 ) =XI*XPR(M) +XR*XPI (M) 

195 

65 

CONTINUE 

196 

GR  = BO 

197 

GI  -0 . 0 

198 

DO  72  11  = 1 , MLS 

199 

GR=GR+B(M)*XPR(M) 

200 

GI=GI+B(M)*XPI(M) 

201 

72 

CONTINUE 

202 

X =ZR*GR“Z I *GI 

203 

GI  =ZR*GI +Z I *GR 

204 

GR=X 

205 

SR=-C1*(GI-FI-FI) 

200 

SI=C1»(GR-FR-FR) 

207 

H2R=GR-SR 

208 

H2I=GI-SI 

209 

GO  TO  317 

210 

120 

FLPS  = 0 

211 

ZM  = DSQRT(ZMSQ) 

212 

ZRT2M  = DSORT(ZM) 

213 

IF  (ZR  .LT.  0 . DO ) GO  TO  125 

214 

ZRT2R  * DSQRT  (0.5D0  • (ZR  * ZM)) 

215 

ZRT2I  « ZI  / (ZRT2R  + ZRT2R) 

216 

Z32R  * ZR»ZRT2R  - ZI»ZRT2I 

217 

Z32I  * ZR»ZRT2I  ♦ ZI*ZRT2R 

218 

GO  TO  130 

219 

125 

ZRT2I  « OSORT  (0 .500  * (ZM  - ZR)) 

220 

IF  (ZI  .LT.  0 . DO ) ZRT2I  » -ZRT2I 

221 

ZRT2R  * ZI  / ( ZR 12 1 ♦ ZRT2I ) 

222 

Z32R  • ZR*ZRT2R  - ZI*ZRT2I 

223 

Z32 I ■ ZR*ZRT2I  ♦ ZI*ZRT2R 

224 

ZM1R  ■ OABS( Z32 I ) 

225 

IF  ( ZM1 R .LT.  TLIM)  GO  TO  130 

226 

R • (TLIM  / ZM1R) 

227 

Z32R  ■ Z32R  • R 
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228 

R - OCBRT(R) 

229 

ZRT2R  • 2RI2R  • R 

230 

ZRT2I  * ZRT2I  * R 

231 

2RT2M  « ZRT2M  * R 

232 

R ■ R * R 

233 

ZM  * ZM  * R 

234 

ZMSO  * ZM* *2 

235 

ZR  « ZR  * R 

236 

ZI  * Z1  * R 

237 

130 

ZRT4R  * DSQRT  (0.5D0  * (ZRT2R  ♦ ZRT2M) ) 

238 

ZRT4I  . O.SOO  • ZRT2I  / ZRT4R 

239 

ZRTM4R  ■ ZRT4R/ZRT2M 

240 

ZRTM4I  > -ZRT4I/ZRT2M 

241 

IF  (ZR  .GT.  0. ) GO  TO  210 

242 

IF  (ZM1R  .IT.  TIIM)  GO  TO  210 

243 

ABK  • ABS(K2) 

244 

IF  (Z321  .GT.  0.)  GO  TO  205 

245 

K1  « K • EXPONT 

246 

K2  « K / EXPONT 

247 

Z32 1 * -TUM 

248 

GO  TO  220 

249 

205 

K3  » K * FXPONT 

250 

HI  ■ K / EXPONT 

251 

Z32I  * TUM 

252 

GO  TO  220 

253 

210 

K2  « C2  ♦ Z32I 

254 

S2  ■ DEXP( K2 ) 

255 

K2  » K»S2 

256 

K 1 « K/S2 

257 

220 

THR  * FPI12  - C2  • Z32R 

258 

STHR  *DS !N( THR ) 

259 

CTHR  .DCOS(THR) 

260 

STP  « -C3*CTHR  +0.5*STHR 

261 

CTP  « -C3*STHR  -0.5*CTHR 

262 

TEMP  • DABS  ( Z32R) 

263 

TEM1  • DABS  ( Z32 I ) 

264 

IF  (TEMP  .LT.  TEM1 ) TEMP  * TEM1 

265 

230 

DO  235  ML  « 1.LA5 

266 

IF  (TEMP  .GT.  ZMLA5(ML) ) GO  TO  250 

267 

235 

CONTINUE 

268 

250 

CONTINUE 

269 

VR  * Z32I 

270 

V I * -Z32R 

271 

CALL  CFR  (YR,  Yl,  F2R.  F2I,  C4 , C5,  ML) 

272 

CPR  • F2R 

273 

CPI  * F2 I 

274 

STORE3=K2* ( ZRTM4R*  F2R-ZRTM4I *F2 1 ) 

275 

ST0RE4*K2«(ZRTM4I«  F2R+ZRTM4R*F2I ) 

276 

H22R  * ST0RE3*C THR- ST0RE4* STHR 

277 

H22I  «ST0RE3*STHR*ST0RE4*CTHR 

278 

IF  (ZR)  280,270,270 

279 

270 

FLQUAO  *0 

280 

GO  TO  300 

281 

280 

IF  (ZI)  290,310,310 

282 

290 

FLQUAD  * 1 

283 

300 

H2R  - H22R 

284 

H2I  > H22I 

285 

GO  TO  317 

286 

3iv 

F LWUMU  * - 1 

287 

YR  « -Z32I 

288 

YI  ■ Z32R 

289 

CALL  CFR  (YR,  YI,  FIR,  FI  I , C4,  C5 

290 

CPR  • FIR 

291 

CPI  * F 1 I 

292 

STORES  = K 1 « ( ZRTM4R*  F 1 R-ZRTM4 1 • F 1 I ) 

293 

ST0RE6*K1*(ZRTM4R«F1 I+ZRTM4I*F1 R) 

294 

H21R=ST0RE5»CTP-ST0RE6«STP 

295 

H21 I*STORE5»STP+STORE6«CTP 

296 

H2R»H2l R*H22R 

297 

H2 I *H21 1+H22I 

298 

317 

IF  (IH  .EQ.  2 )GO  TO  80 

299 

60 

IF  (FLPS  .NE.  1 ) GO  TO  320 

300 

70 

HI R * GR+SR 

301 

HI  I • GI+SI 

302 

GO  TO  362 

303 

320 

IF  (FLQUAD  . LT.  0 ) GO  TO  340 

304 

330 

YR  • -Z32I 

305 

YI  ■ Z32R 

306 

CALL  CFR  (YR,  YI,  FIR,  FI  1 , C4,  C5 

307 

340 

ST0RE7  = K1» ( ZR  TM4R • F 1 R-ZRTM4 1 *F 1 I ) 

308 

ST0RE8  = K 1 * ( ZRTM4 1 * F 1 R-f  ZRTM4R* F 1 I ) 

309 

HI  1 R=STORE7*C THR+STORE8*STHR 

310 

HI  1 I*STORE7*(-STHR)+STORE8»CTHR 

311 

IF  ( FLQUAD  .LE.  0)  GO  TO  360 

312 

350 

ST0RE9=K2» ( ZR  TM4R*  F2R~ ZRTM4 1 *F2 1 ) 

313 

STOR10=K2*(ZRTM4I*F2R+ZRTM4R*F2I) 

314 

H12R  » STORE9*CTP+ST0R10*STP 

315 

HI  2 1 * STORE9*(~STP)+STOR10*CTP 

316 

H1R  = HI 1R+H12R 

317 

HI  I * HI  1 1+H12I 

318 

GO  TO  362 

319 

360 

HI R * HI 1R 

320 

HI  I * Hill 

321 

362 

IF  ( IH  . EQ.  1 ) GO  TO  999 

322 

80 

IF  (FLPS  .NE.  1)  GO  TO  380 

323 

90 

FPR  * CO 

324 

FPI  * 0.0 

325 

00  9?  M * 1 , MLS 

326 

FPR*FPR+C(M)»XPR(M) 

327 

92 

FPI«FPI+C(M)*XPI(M) 

328 

X «-( STORE 1*FPR-ST0RE2*FPI) 

329 

FPI*-(ST0RE1*FPI4-ST0RE2*FPR) 

330 

FPR  ■ X 

331 

GPR  » DO 

332 

GPI  * 0.0 

333 

00  94  M * 1 , MLS 

334 

GPR«GPR+D(M) *XPR(M ) 

335 

94 

GPI «GPI +D(M)*XPI(M) 

336 

SPR--C1 «(GPI-FPI-FPI) 

337 

SPI»C1»(GPR-FPR-FPR) 

33(i 

H2PR*GPR“SPR 

'J3  9 

H2PI-GPI-SPI 

340 

GO  TO  414 

341 

380 

YR  ■ Z32I 
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342 

Y I • -Z32R 

343 

CALL  CFR  (YR,  YI,  G2R , G2I,  D4,  D5,  ML) 

344 

STOB11  » K2» ( ZRT4R*G2R_ZRT4I *G2I ) 

345 

STOR12  « K2« ( ZRT4R*G2I+ZRT4l *G2R ) 

346 

H22PR«STOR1 1 •STHR+ST0R12*CTHR 

347 

H22P1«ST0R11* (-CTHR)  •fSTOR12*STHR 

348 

390 

IF  (FLOUAD  .LT.  0)  GO  TO  410 

349 

400 

H2PR  ■ H22PR 

350 

H2PI  * H22PI 

351 

GO  TO  414 

352 

410 

YR  ■ -Z32I 

353 

YI  ■ Z32R 

354 

CALL  CFR  (YR,  YI,  G1R,  Gil,  D4,  05,  ML) 

355 

ST0R13  > K1*(ZRT4R*G1 R-ZRT4I »G1 I ) 

356 

STORM  • K1*(ZRT4R»G11+ZRT4I»G1R) 

357 

H21PR»STOR13«(-STP)  -ST0R14»CTP 

358 

H21PI«ST0R14*(-STP)  +ST0R1 3*CTP 

359 

H2PR  > H2IPR+H22PR 

360 

H2PI  • H21P1+H22PI 

361 

414 

IF  ( IH  .EQ.  2)  GO  TO  999 

362 

100 

IF  (cLrS  .NE.  1)  GO  TO  420 

363 

110 

HI  PR  « GPR+SPR 

364 

HI  PI  « GPUSPI 

365 

GO  TO  999 

366 

420 

IF  ( FLQUAD  .LT.  0)  GO  TO  440 

367 

430 

YR  « -Z32I 

368 

YI  • Z32R 

369 

CALL  CFR  (YR,  YI,  G1R,  G1 I , D4 , 05.  ML) 

370 

440 

ST0R15  * K1 * ( ZRT4R»G1 R -ZRT4I*G1I) 

371 

ST0R16  » HI  * ( ZRT4R*G1 1 +ZRT4I*G1R) 

372 

H 1 1 PR  . ST0R1 5*STHR  -ST0R16*CTHR 

373 

HI  1 PI  » ST0R15«CTHR  +STOR 1 6*STHR 

374 

450 

IF  (FLOUAD  .GT.  0)  GO  TO  470 

375 

460 

HI  PR  « H1 1 PR 

376 

H 1 P I • H1 1 PI 

377 

GO  TO  999 

378 

470 

ST0R17  « K2*  ( ZRT4R»G2R  -ZRT4I*G2I) 

379 

STOR 1 8 > K2*  ( ZRT4R  »G2 1 •fZRT4 1 »G2R ) 

380 

H12PR  « ST0R1 7* (-STP)  ♦ST0R18»CTP 

381 

H12PI  « ST0R1 7*(-CTP)  -ST0R18*STP 

382 

HI  PR  » H12PR+H1 1 PR 

383 

HI  PI  ■ H12PI+H11PI 

384 

999 

CONTINUE 

385 

RETURN 

386 

END 

•PRT.S 
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1 SUBROUTINE  CFR(X,  V,  SR,  SI,  A,  B,  M) 

2 IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

3 DIMENSION  Ail).  B( 1 ) 

4 SR  • O.DO 

5 SI  ■ O.DO 

6 DO  tO  d ■ t.M 

7 1 • M - J ♦ 1 

8 TEMR  « X ♦ SR  ♦ B(  I ) 

9 TEMI  » Y + SI 

10  TEMP  • A(I)  / (TEMR**2  ♦ TEMI»*2) 

11  SR  ■ TEMR  * TEMP 

12  SI  ■ -TEMI  • TEMP 

13  10  CONTINUE 

14  SR  ■ SR  4-  1 .DO 

15  RETURN 

16  END 


t 

SUBROUTINE  HZERO(PARTR, PARTI ) 

2 

IMPLICIT  DOUBLE  PRECISION( A-H.O-Z) 

3 

DOUBLE  PRECISION  PARTI,  PARTR,  K2,  IMZ12 

4 

COMMON  /AHZERO/  HZEROR ,H2ER0I 

S 

02  ■ PAR TR* • 2 ♦ PART  I a»2 

• 

8 K2  • . 7978845608* EXP (PART  I ) 

7 

D ■ SORT (D2) 

8 

RLZI2  • (SORT( (PARTR  ♦ D)/2.))/D 

9 

I F ( D - PARTR)  9,9,10 

• 

10 

9 IMZ12  > 0 

11 

GO  TO  11 

12 

10  IMZ12  « ( SORT ( (D  - PARTR)/2. ) )/D 

13 

11  COST  ■ COS(PARTR) 

14 

SINT  a-SIN(PARTR) 

15 

HZEROR  ■ K2* ( RLZ1 2*C0ST  - IMZ12*SINT) 

18 

HZEROI  > K2* ( IMZl 2* COST  + RLZ12*SINT) 

17 

RETURN 

18 

END 
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APPENDIX  B:  SAMPLE  RUN 


This  appendix  gives  a brief  discussion  of  the  input-output,  then  lists  an  input  deck 
and  shows  the  resulting  output.  The  input  deck  is  really  three  separate  runs  that  are  stacked 
to  run  consecutively.  The  input  to  a single  run  consists  of  several  parts  given  in  table  Bl. 
The  table  gives  the  number  of  cards  and  the  location  of  the  FORTRAN  input  statements  in 
Program  MAIN.  The  last  three  of  these  are  open-ended.  That  is,  more  modes,  receiver 
depths,  or  ranges  are  read  in  until  a blank  card  specifies  the  end  of  the  list.  A blank  range 
card  sends  the  program  back  to  the  beginning.  A negative  range  sends  the  program  back  to 
read  a new  source  and  new  receivers  after  reading  another  key  card.  The  program  halts 
when  a blank  “n  and  freq”  card  is  encountered. 

Table  B2  gives  most  of  the  functions  of  the  key  card  by  which  integers  are  read  into 
control  keys  1-9.  Some  of  these  will  require  additional  information,  which  is  read  in  imme- 
diately following  the  key  card. 

The  output  of  the  program  is  usually  printed  through  FORTRAN  print  statements. 
Cards  are  also  punched  when  key  5 is  10  or  key  8 is  2,  3,  or  4.  In  the  first  case  each  card 
contains  a complete  eigenvalue  that  can  be  read  into  future  runs. 

When  key  8 = 2,  propagation  losses  for  25  consecutive  ranges  per  card  are  punched 
for  each  receiver  depth,  with  a maximum  of  5 receiver  depths.  The  losses  can  be  read  into  a 
plot  program  with  a format  of  (5X.25F3.1).  Each  loss  must  then  be  subtracted  from  140. 
This  format  allows  losses  to  tenths  of  a dB  from  40.1  to  140.0  dB. 

When  key  8 is  3,  losses  for  up  to  26  receiver  depths  are  punched  on  one  card  for 
each  range.  These  cards  can  be  used  in  a contour  plotting  program.  They  can  be  read  with  a 
format  of  (26F3.1 ) and  must  also  be  subtracted  from  140. 


Table  Bl.  Input  cards  to  the  normal  mode  program. 


Input 

Function 

Number  of  Cards 

Location 
in  Program 
MAIN 

Control  keys 

Selects  options 

1 or  more 

37-65 

n and  freq 

Determines  number  of  layers  and  fre- 
quency; also  halts  program 

1 

66 

Profile 

Specifies  depths,  sound  speeds,  gradients, 
attenuations,  and  densities 

7 

71-85 

Modes 

Searches  for  or  specifies  modes 

1 or  more 

224 

Source  and  receivers 

Specifies  a source  depth  and  one  or 
more  receiver  depths 

2 or  more 

461-463 

Ranges 

Specifies  a sequence  of  ranges;  also 
directs  continuation 

1 or  more 

616 
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Table  B2.  Functions  of  the  control  keys. 


Key 

Setting 

Effect 

Function  Affected 

1 

0 

No  output 

Print 

Depth  functions 

Plot  on  printer 

2 

0 

Print  losses 

Propagation  losses 

1 

Contour  on  printer 

3 

0 

No  output 

Reflection  coefficients 

1 

Print  all  interfaces 

k > l 

Print  interface  k 

4 

k >0 

Change  levels  and  symbols 

Contour  on  printer 

5 

0 

Sum  only  those  given 

Number  of  modes 

1 

Add  to  existing  sum 

k + 10 

Punch  modes  on  cards 

6 

0 

Long  print 

Steps  in  mode  iteration 

1 

Short  print 

2 

Shortest  print 

7 

0 

Phased  addition 

Mode  sum 

1 

Random-phase  addition 

8 

1 

Change  T-lim 

2 

Punch  losses  for  up  to  S 

Loss  plot  input 

receivers 

3 

Punch  losses  for  up  to  26 

Contour  plot  input 

receivers 

9 

0 

No  effect 

Number  of  modes 

k 

Use  only  1st  k modes 

The  first  profile  in  the  input-output  is  a surface  duct.  100  m deep.  For  the  500  Hr 
frequency,  3 modes  are  found  by  searching  from  a phase  velocity  of  1 520.5-1 523  m/s.  Two 
additional  modes  are  found  by  extrapolation.  Forty  receiver  depths  are  then  specified  from 
3 to  120  m,  and  propagation  loss  contours  are  drawn  for  a source  at  a depth  of  20  m.  The 
modes  are  added  in  random  phase,  and  loss  contours  of  80. 90.  and  100  dB  are  requested  to 
be  represented  by  the  symbols  8,9,  and  0.  A negative  range  then  causes  the  program  to  read 
new  control  keys,  source  and  receivers.  The  depth  functions  are  then  printed  out  as  ampli- 
tudes and  phase  angles  and  propagation  losses  are  computed. 

The  second  profile  consists  of  two  negative  gradients  over  two  layers  of  sediments  in 
shallow  water.  A velocity  discontinuity  exists  at  the  top  of  each  sediment  layer.  Negative 
numbers  in  the  input  for  the  attenuation  at  the  bottom  of  the  sediment  layers  serve  as  flags 
to  request  that  the  gradients  at  the  top  of  the  layers  have  no  imaginary  parts  and  that  the 
attenuation  at  the  bottom  will  be  whatever  results  from  this.  The  change  in  ImC  from  37.9 
to  23.7  in  the  deeper  sediment  layer  indicates  that  the  attenuation  changed  by  about  60 
percent  through  the  layer. 


A final  layer  of  negative  sound  speed  gradient  must  always  be  added.  A gradient  of 
-0. 1 is  chosen  here  for  the  top  of  this  layer. 

The  first  three  modes  are  determined  by  reading  in  approximate  values.  The  magni- 
tudes of  the  depth  functions  are  plotted  on  a log  scale  at  2-m  depths  from  30  to  80  m. 
Reflection  coefficients  are  computed  at  interface  2,  which  is  the  water-sediment  interface. 

The  final  profile  is  a deep-water  profile  with  a 40-m  deep  sediment  layer.  The  atten- 
uation increases  from  2 dB/km  to  2.5  dB/km  through  this  layer.  The  first  mode,  the  first 
bottom-reflected  mode,  and  a higher  bottom-reflected  mode  are  found.  Each  step  of  the 
mode  iteration  is  printed  out.  Reflection  coefficients  are  again  computed.  The  amplitudes 
and  phase  of  the  depth  functions  are  printed  out  at  500-m  depth  and  at  each  even  1 000-m 
depth  for  a 100-m  source  depth. 

On  the  last  two  profiles,  a much  larger  set  of  modes  is  required  to  compute  correct 
propagation  losses. 

The  sample  run  given  here  required  6 seconds  on  a UNIVAC  1110,  Exec  8 operating 
system.  The  cost  of  the  run  was  $1 .20. 


L 
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NORMAL*MODE(0) .INPUT 


t 

2 

3 

4 

5 

e 

7 

8 
9 

10 

11 

13 

13 

14 


50 

51 
53 

53 

54 

55 

56 


INPUT  DECK  STARTS  AT  LINE  3,  ENOS  AT  LINE  65. 

123456789  123456789  123456789  123456789  123456789  123456789  123456789  123458789 


100. 

0 

098 

2 

1520. 

0 

.017 

0 

0 

0 


500. 


90. 


100. 


-.1 


1 

80. 


2 

-1000. 


15 

-1520.5 

1523. 

16 

-1  . 

17 

0 

18 

10. 

19 

3. 

120. 

3. 

20 

0 

21 

4000. 

100000. 

4000. 

22 

-1  . 

23 

1 

24 

10. 

25 

30. 

120. 

30. 

26 

0 

27 

5000. 

100000. 

5000. 

28 

0 

29 

2 

2 

1 

30 

5 1500. 

31 

0. 

51  . 

73. 

73.3 

373.3 

32 

1542.2 

1536.8 

1606.45 

1684. 

33 

1523.42 

34 

1 .5 

1.5 

-.1 

35 

.12 

.73 

.73 

36 

-1 . 

-1 . 

37 

1.68 

1 .91 

1.91 

38 

1527.18 

.16 

39 

1530.64 

.13 

40 

1533.49 

.11 

41 

0 

42 

60. 

43 

30. 

80. 

2. 

44 

0 

45 

0 

46 

1 

6 

47 

8 100. 

48 

55. 

146. 

402. 

960. 

49 

1544.9 

1542.6 

1517.9 

1495.0 

1483.: 

30 

2 


2286. 

1497.8 

1533.4 


4390. 
1541 ,7 

1 . 

.02 
.025 
1 .54 


-1483.5 

-1533.4 


.1 


1484.5 

1534.4 


10 

10 


4430. 


-.  1 

.025 

2.5 
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APPENDIX  C:  HANKEL  FUNCTION  PARAMETERS 

This  appendix  gives  the  FORTRAN  statements  for  two  programs  associated  with  the 
modified  Hankel  functions.  Program  PWRTRN  computes  the  power  series  coefficients,  dm, 
from  eq  (57),  then  determines  the  truncation  points  from  eq  (59).  The  truncation  points 
for  the  other  three  sets  of  coefficients  can  be  determined  by  changing  line  9.  Different  com- 
puter word  lengths  can  be  accommodated  by  changing  line  16. 

The  second  program,  CFC,  determines  the  asymptotic  series  coefficients  Cm  from 
eq  (72),  then  determines  the  continued  fraction  coefficients  as  indicated  by  eq  (81M83). 
The  second  set  of  coefficients  can  be  determined  by  changing  the  4 in  line  1 1 to  a 1 6. 
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PROGRAM  PWRTRN 

C *•  THIS  PROGRAM  DETERMINS  TRUNCATION  POINTS  FOR  THE  POWER  SERIES. 
IMPLICIT  DOUBLE  PRECISION A-H.O-Z) 

DIMENSION  0(50),  ALOGD(50) 

0(1)  . 1. 

A LOGO) 1 ) a 0. 

P a 3 , 

OQ  50  I a 2,50 

0(1)  - 0(1-1)  / P / (P-2) 

P a p «•  3. 

IF  (0(1)  . LE  0.)  GO  TO  50 
ALOGO(I ) a AL0G1 0( 0(1)) 

50  CONTINUE 

PRINT  60,  D,  ALOGO 
60  FORMAT  (10E12.6) 

OH  a 18. 

M a 1 

DO  10  K a 2.50 
30  P a M - K 

2 a (ALOGD(K)  - ALOGD(M)  ♦ DH)  / 3.  / P 
IF  (P  .GT.  -1.1)  GO  TO  20 
A • ALOGD(M)  - ALOGO(M-M)  - 3.  a Z 
IF  (A  .GT.  0. ) GO  TO  20 
M • M ♦ 1 
GO  TO  30 
20  L ■ K - 1 
MM  a M - 1 

A Z a EXP  (Z  • 2.3025851) 

AZSO  » AZ  * AZ 

PRINT  40,  L,  MM,  Z.  AZ,  AZSQ 
40  FORMAT  (215,  4E15.8) 

10  CONTINUE 
ENO 


90 


PROGRAM  CFC 

C **  THIS  PROGRAM  COMPUTES  A SET  OF  SERIES  COEFFICIENTS  AND  THEN 
C ••  COMPUTES  THE  CORRESPONDING  CONTINUED  FRACTION  COEFFICIENTS. 
IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

DIMENSION  COEF(21 ,23,3).  CHECK(20),  C(82),  S(10).  A(20).  B(20) 
C(1)  - 1. 

BOTTOM  ■ 1. 

TOP  ■ 1 . 

DO  2 I • 1 .45 
X . 48  * I 

V » 9 * (I  + I - 1 )*»2  - 4 
C(I+1)  * C(I)  * V / X 
2 CONTINUE 

PRINT  20,  (C(I),  I ■ 1.40) 


PRINT  20,  (C(I),  I ■ 1.40) 
FORMAT  (5G20.9) 

FORMAT  (/) 

DO  100  I • 1,11 
COEF( 1,1,3)  » 0. 

COEF ( I , 1+1 ,3)  » 0. 

COEF ( I , 1+2 ,3 ) » 0. 

CONTINUE 
A ( 1 ) • C(2) 

COEF( 2 ,2,3)  * 1. 

DO  140  I « 3,21 

DO  110  J « 2.1 

COEF(I.J.I)  ■ COEF ( 1-1 , J , 3) 

COEF ( I , J ,2 ) « COEF ( 1—2 , J • 3 ) 

COEF ( I , J , 3 ) - COEF( 1-1 .J-1 ,3! 

CONTINUE 

IF  (I  .ED.  3)  GO  TO  150 
CON  « 0. 


AT  » 0. 

BT  m 0. 

K « I - 3 
DO  120  J ■ 3,1 
K ■ K + 1 
CON  » C ( K ) • C 


CON  » C ( K ) • COEF( I.J-1,3)  + CON 
AT  « C ( K ) * COEF(I , J-1 ,2)  ♦ AT 
BT  « C ( K ) • COEF(I , J-1 ,1)  + BT 
CONTINUE 

PRINT  160,  CON,  AT,  BT 

CHECK( 1-2)  « BT 

A ( 1-2 ) * -(CON  + C ( K+1 ) ) / AT 

CONTINUE 

CON  « 0. 

AT  ■ 0. 

BT  » 0. 

K » I - 2 
DO  130  J « 3.1 
K « K + 1 

CON  « C ( K ) • COEF( I , J-1 ,3)  + CON 
AT  ■ C(K)  . COEF ( I , J-1 ,2)  ♦ AT 
BT  « C(K)  • COEF ( I , J-1 , 1 ) + BT 
CONTINUE 

PRINT  160,  CON,  AT,  BT 
PRINT  11 


E 


■j 


* 
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APPENDIX  D:  MODE  FOLLOWER  PROGRAM  IN  FORTRAN 


The  FORTRAN  statements  of  the  Mode  Follower  program  are  given  here.  This  is 
the  main  body  of  the  program.  The  following  auxiliary  subroutines  from  appendix  A are 
required:  SETUP,  DET,  HANKEL,  and  CFR. 
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AD-A072  201  NAVAL  OCEAN  SYSTEMS  CENTER  SAN  DIE60  CA  F/G  17/1 

UNDERWATER  SOUNO  PROPAgATION-LOSS  PROGRAM.  COMPUTATION  BY  NORMA— ETC (U) 
MAY  79  D F GORDON 

UNCLASSIFIED  NOSC-TR-393  NL 


2 » 2 

AO 

4072.  Of 


1 

a 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 
58 


PROGRAM  MFOLLO 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  /INPUT/  Z(10) ,N, OMEGA, V, VI ,GCU( 10) . GSO( 10) ,CAY( 1 0) , LAMBDA, L 
1 AM80I , G( 1 0) , RHO( 1 0 ) ,GI ( 1 0 ) , GSQI (10), CAY  1(10) 

COMMON  /DETMNT/  A( 21 ,4) ,0( 21 ,4) 

REAL  INCA,  INCB,  INCC,  INCD,  INCE,  LAMBDA,  LAMBDI 
DIMENSION  T( 4 ) , PV(4),  W(8),  WI ( 8) , CB(10),  CBI(IO),  C( 10) . 

1 CAY  S0(10).  GAMMA ( 10)  , DPK(IO),  GCUI(IO),  CI(10),  CR( 1 0) ,PVI (4) 

2 , CAYSQI ( 10 ) , SR( 4 ) , SI(4) 

CHNG  * 1 . / 8192. 

CHNGI  * 0. 

4 CONTINUE 

Cii  1(8  - TOTAL  STEP  LIMIT,  K1  , K2  PRINT  KEYS,  K3  * i ivEEPS  SamE  PHuriLC 


C** 


10 


30 


FOR  NEXT  RUN. 

READ  10,  KO,  K1 , K2.  K6,  K3,  TLIM,  BLIM,  RATIO,  EX 
PRINT  10, KO,  K1 , K2 , K6,  K3,  TLIM,  BLIM,  RATIO,  EX 
FORMAT  (514,  4E10.1) 

IF  (TLIM  . EQ . 0.)  TLIM  • 1.E-5 
IF  (BLIM  .EQ.  0.)  BLIM  * 1 . E-2 
IF  (EX  .EQ.  0. ) EX  - 28. 

RLIM  * 1 0 . **EX 

IF  (RATIO  .EQ.  0.)  RATIO  « 2. 

IF  (KO  .EO.  0)  KO  * 300 
IF  (K3  .NE.  0)  GO  TO  128 
READ  1240,  N.FREQ 


.ATTEN 

C**  STOP  IF  N = 0.  THIS  IS  THE  ONLY  PROGRAMED  STOP. 

IF  (N.EQ.O)  GO  TO  1200 
PRINT  1250,  N.FREQ 

C**  PARAMETERS  READ  IN  BELOW  ARE  THOSE  AT  THE  TOP  OF  EACH  LAYER. 
C**  REAO  IN  VELOCITIES. 

READ  1260,  (C( I ) , I *1 ,N) 

PRINT  1280.  (C(I).I-I.N) 

C**  READ  IN  DEPTHS. 

READ  1260,  ( Z ( I ) , I * 1 ,N) 

PRINT  1280,  (Z(I),I*1,N) 

C**  READ  IM  GRADIENTS 

READ  1260,  ( GAMMA ( I ) , I • 1 , N ) 

PRINT  1280,  (GAMMA ( I ) , 1*1 ,N) 

C*  * READ  IN  ATTENUATION  FACTOR  IN  LOSS  PER  KILOMETER. 

READ  1260,  (DPK(I) ,1*1 ,N) 

PRINT  1280,  ( DPK( I ) , I *1 ,N) 

C**  READ  IN  DENSITIES  (BLANK  INPUT  IMPLIES  SEA  WATER  DENSITY). 
READ  1260,  (RHO( I) , 1*1 ,N) 

PRINT  1280,  (RHO( I ) , 1-1 ,N) 

CONTINUE 
NUMBER  * 1 
JX  - 0 

NX  * VARIABLE,  NY  * LAYER  NUMBER,  NZ 
READ  119,  NX,  NY,  NZ,  PK,  VALL.OP,  V,  VI,  STEP,  STEPI 
PRINT  21,  NX,  NY,  NZ,  PK,  VALL.DP,  V,  VI,  STEP,  STEPI 
FORMAT  i 312 , 4X,  7010.2) 

FORMAT  ( 1 0H  VARIABLE  .12,  10H  LAYER  NO 
* 12,  / 7G15.5) 

PK  - PK  - DP 

C**  START  NEW  CYCLE  BY  INCREMENTING  VARIABLE. 

107  PK  ■ PK  ♦ DP 


128 


C** 


119 

21 


CONTINUITY 


I2.12H  CONT!NHtrv 


94 


57 

IF  (OP)  108,999,109 

88 

C**  ( 

CHECK  IF  OESIRED  LIMIT  OF  VARIABLE  HAS  BEEN  REACHED. 

89 

108 

IF  (PK  . LT . VALL)  GO  TO  3 

60 

GO  TO  133 

ei 

109 

IF  (PK  .GT.  VALL)  GO  TO  3 

62 

133 

GO  TO  (131, 101, 102, 103, 104, 105), NX 

63 

131 

FREO  ■ PK 

64 

GO  TO  106 

68 

101 

C(NY)  ■ PK 

66 

IF  (NZ  .NE.  0)  GO  TO  106 

67 

134 

IF  (NY  .EQ.  N)  GO  TO  135 

68 

GAMMA(NY)  • 0. 

69 

IF  (NY  .LT.  2)  GO  TO  106 

70 

138 

GAMMA ( NY- 1 ) ■ 0. 

71 

GO  TO  106 

72 

102 

Z(NY ) ■ PK 

73 

IF  (NZ  .EQ.  1)  GO  TO  106 

74 

IF  (NY  .LT.  N)  GO  TO  134 

78 

IF  (NUMBER  .EQ.  1)  GO  TO  106 

76 

C ( NY ) - 0. 

77 

GO  TO  106 

76 

103 

GAMMA(NY)  « PK 

79 

IF  (NZ  . NE.O)  GO  TO  106 

80 

J • NY  ♦ 1 

61 

DO  121  I ■ J.N 

82 

C(I)  - 0. 

83 

121 

CONTINUE 

84 

104 

DPK(NY)  • PK 

88 

GO  TO  106 

86 

105 

RHO(NY)  « PK 

87 

106 

CONTINUE 

68 

89 

C*»  i 

COMPLETE  PROFILE  *• 

90 

DO  100  1-1 , N 

91 

C**  ! 

SET  UNSPECIFIED  DENSITIES  TO  1.02  (SEA  WATER). 

92 

IF  (RHO(I).NE.O. ) GO  TO  40 

93 

RHO(I )«1 .02 

94 

40 

IF  (I .EQ.1 ) GO  TO  50 

95 

C** 

COMPUTE  VELOCITY  AT  BOTTOM  OF  PREVIOUS  LAYER. 

96 

TEMP* Cl (1-1 ) **2 

97 

TEMDR-C( 1-1  )--2 

96 

TEMDI * ( TEMDR+TEMDR+TEMDR-TEMP ) *CI ( 1-1 ) 

99 

TEMDR* ( TEMDR-TEMP-TEMP-TEMP) *C( 1-1 ) 

100 

TEMP* ( GAMMA ( 1-1 ) +GAMMA( 1-1 ) ) * ( Z( I )-Z( 1-1 ) )-C( 1-1 ) 

101 

TEMDEN-TEMP* -2+C I ( 1-1 )*»2 

102 

TEMI* ( TEMDI *CI( 1-1 ) -TEMDR* TEMP )/TEMDEN 

103 

TEM1 I * (-TEMDI-TEMP-TEMDR-CI ( 1-1 ) )/TEMDEN 

104 

CB( I )*SQRT ( . 5* ( TEMI  •fSQRT( T EMI **2+TEMl 1**2) ) ) 

105 

CBI ( I ) -TEMI I/(C8( I)+CB( I )) 

106 

SO 

IF  (C(I).NE.O)  GO  TO  60 

107 

c** 

IF  VELOCITY  WAS  UNSPECIFIED  USE  VELOCITY  AT  BOTTOM  OF  PREVIOUS  LAYER 

108 

C( I )»CB( I ) 

109 

60 

IF  (OPKU1.NE.0.  ) GO  TO  70 

110 

CI(I)-0. 

111 

GO  TO  60 

112 

C** 

IF  ATTENUATION  IS  TO  BE  APPLIED  TO  A LAYER,  COMPUTE  COMPLEX  VELOCITY 

113 

C*«  l 

KEEP  ABSOLUTE  C EQUAL  TO  GIVEN  REAL  C FOR  SIMPLICITY. 

95 


I 

i 

| 


114 

70 

TEMP«S4574T-*FREQ 

118 

TEMOI *DPK( I ) »C( I ) 

lie 

TEMDR.TEMP**2+TEMDI**2 

117 

Cl(l)«TEMDI«TEMP*C( I )/TEMOR 

lit 

C( I )«TEMP**2«C( I ) / T EMOR 

lit 

60 

IF  ( GAMMA ( I ) .Nt . 0 . ) GO  TO  100 

120 

IF  ( 1 . EQ.N ) GO  TO  90 

121 

C** 

COMPUTE  GRADIENT  IF  NOT  GIVEN. 

122 

GAMMA ( I ) . ( C ( I'M  )**2-C(I )«*2)*C(I)/(2.*C( 1*1 )**2« (Z( 1+1 )— Z( I ) ) ) 

123 

IF  (l.EQ.N)  GO  TO  90 

124 

GO  TO  100 

128 

C*« 

REDUCE  LAYERS  6Y  ONE  IF  FINAL  POINT  ONLY  DEFINES  GRADIENT  IN  LAST  LAVER. 

126 

90 

N«N-1 

127 

100 

CONTINUE 

126 

129 

C*» 

COMPUTE  USEFULL  0UANT1ES  •* 

130 

OMEGA. 6 . 283185307’ FREO 

131 

DO  120  I . 1 ,N 

132 

TEMP.C( I )**2+CI(  ! )**2 

133 

CAY( II .OMEGA«C( I )/TEMP 

134 

CAYI( I).-OMEGA»CI( D/TEMP 

138 

CAYSQ( I ) «C A Y ( I ) * .2-CAYI ( I ) **2 

136 

CAYSOU I ).2.*CAY(I )*CAVI( I ) 

137 

TEMDR.-2. *GAMMA( I )*CAYSQ(I) 

138 

TEMDI.-2. .GAMMA ( I ) *CAYSQI ( 1 ) 

139 

GCU(  I ) « ( TEMDR*C(  I )-M£MD  I*CI(I))/TEMP 

140 

GCU1( I ).(TEMDI*C( I )-TEMOR*CI( I ))/TEMP 

141 

TEMP«EXP( ALOG( GCU( I ) * »2+GCUl ( I ) **2 )/6 . ) 

142 

GI(I).TEMP>SIN(ATAN2(GCUl(l),ABS(GCU(I)))/3. ) 

143 

G(I)*SQRT(TEMP»*2-GI(I)«*2) 

144 

IF  ( GAMMA ( I ) . LT . 0 . ) GO  TO  110 

148 

G(  I ) **G( I ) 

146 

110 

GI(I).-GMI) 

147 

£ * * 

.*.»  IS  LAVS'?  STRENGH  PARAMETER  USED  ONLY  TO  COMPARE  WITH  OTiiER  MODS 

148 

XMI.-GK I )*(2( 1+1 )-2(I) ) 

149 

XM«-G( T)*(2(I+1)”E(I)) 

ISO 

GSOI (1)>2.*G(I)*GI(I) 

181 

120 

GSO( I )»G( I ) **2-GI ( I )**2 

182 

IF  (JX  .GT.  0)  GO  TO  113 

153 

C** 

GO  TO  INITIAL  3 STEPS  OR  TO  THE  STANOARD  STEP. 

154 

IF  (NUMBER  - 4)  71,111,122 

155 

71 

CALL  SETUP 

156 

CALL  DETNT(N,DET,DETI) 

157 

VEL.DET 

158 

VELI.DETI 

159 

DELTA. STEP 

160 

DELT1.STEPI 

161 

IF(OELTA.NE.O. )G0  TO  250 

162 

1F(DELTI.E0.0.)DELTA«.01 

163 

280 

SIZE2. 100. 

164 

IF  (K6.LT. 3)  PRINT  1320,  V , VI , DET ,DETI , A ( 21 ,4) ,Q(21 ,4) 

168 

C** 

ITERATE  FOR  MOOE  UP  TO  7 STEPS. 

166 

DO  310  J.1 ,12 

167 

V.V+DELTA 

168 

VI.VIfOELTI 

169 

C** 

00  NOT  PERMIT  IMAGINARY  PART  TO  BECOME  NEGATIVE. 

170 

IF  (VI)  260,270,260 

96 
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171 

260 

DELTI-OELTI-VI 

i7a 

270 

Vial. E-IB 

173 

C** 

SET  UP  DETERMINANT  FOR  PHASE  VELOCITY  V 

♦ VI 

174 

280 

CALL  SETUP 

175 

CALL  DETNT  (N.DET.DETl) 

178 

IF  (K6.NE.1)  GO  TO  300 

177 

PRINT  1330.  V.  VI,  GET,  DETI,  SLR,  SLI 

178 

300 

TEMNR  - DELTA 

179 

TEMNI  > DELT1 

180 

TEMDR-VEL-OET 

181 

TEMDI-VELI-DETI 

182 

TEMDEN*TEM0R*TEM0R>TEMD1*TEM0I 

183 

IF  (TEMDEN.EO.O. ) GO  TO  320 

184 

TC;.l«;iw-7Cr«NR*TEM0R+TEMNI  *TEMDI 

185 

T EMI NU» TEMNI •TEMDR-TEMNRaTEMDI 

186 

SLR  ■ TEMRNU  / TEMDEN 

187 

SLI  a TEMINU  / TEMDEN 

188 

IF  (J  .GT.  3)  GO  TO  125 

189 

• 

SR(4-NUMBER)  ■ SLR 

190 

SI (4-NUMBER)  - SLI 

191 

125 

DELTA  ■ DET  a SLR  - DETI  • SLI 

192 

DELTI  . DET  • SLI  ♦ DETI  • SLR 

193 

size*delta*o£lta+delti*delti 

194 

C*« 

DISCONTINUE  ITERATION  AFTER  2ND  STEP  IF 

CORRECTION 

195 

c** 

PREVIOUS  STEP. 

19C 

IF  ( ( SIZE . GT . S IZE2) .AND. ( J .GT .2) ) 

GO  TO  320 

197 

SIZE2-SIZE*2. 

198 

VEL«DET 

199 

VEL I «DET I 

200 

310 

CONTINUE 

201 

320 

CONTINUE 

202 

51 

PV( 4-NUMBER)  • V 

203 

PV I (4-NUMBER) ■ VI 

204 

NUMDER  a NUMBER  ♦ 1 

205 

GO  TO  107 

206 

C** 

START  STANDARD  STEP,  EXTRAPOLATE  PHASE 

VELOCITY  AND 

207 

111 

INCA  a DP 

208 

INCB  a DP 

209 

INCC  a OP 

210 

122 

INCD  a -INCB  - INCC 

211 

T ( 1 ) a INCB  • INCD 

212 

T( 2)  a INCB  * INCC 

213 

T ( 3 ) - INCD  • INCC 

214 

DO  112  IS  ■ 1,3 

215 

N( IS+4)  - -SR(IS)  ! TC IS) 

216 

MI ( IS  ♦ 4)  • -SI(IS)  / T( IS) 

217 

W(  IS)  a -PV(IS)  / T ( IS) 

218 

112 

MI ( IS)  . -PVI(IS)  / T( IS) 

219 

113 

INCD  a INCA  ♦ INCB 

220 

INCE  a INCD  ♦ INCC 

221 

T ( 1 ) a INCD  a INCE 

222 

T(2)  a INCA  a INCE 

223 

TO)  a INCA  a INCD 

224 

SLOP  a 0. 

225 

SLOPI  a 0. 

226 

SUM  a 0. 

227 

SUMI  a 0. 

97 


238 

00  114  IS  ■ 1,3 

229 

SLOP  - SLOP  ♦ W( IS  44)  • T( IS) 

230 

SLOPI  ■ SLOPI  ♦ MI  ( ISM)  • T(  IS) 

231 

SUM  ■ SUM  ♦ M( IS)  * T( IS) 

232 

114  SUMI  ■ SUMI  ♦ Mil  IS)  • T( IS) 

233 

V ■ SUM 

234 

VI  ■ SUMI 

235 

CALL  SETUP 

236 

CALL  OETNT  (N.DET.DETI) 

237 

C**  EVALUATE  DETERMINANT  AT  THE  EXTRAPOLATED 

POINT. 

238 

VEL  - GET 

239 

VELI  - DETI 

240 

C«*  ITERATE  FOR  THE  ROOT  USING  EXTRAPOLATED  SLOPE. 

241 

DELTA  ■ DET  • SLOP  - DETI  • SLOPI 

242 

DELTI  ■ DET  • SLOPI  ♦ DETI  • SLOP 

243 

IF  (K1  .EQ.  1)  PRINT  1330,  V,  VI,  DET. 

DETI. 

DELTA,  DELTI 

244 

V • V ♦ DELTA 

245 

VI  ■ VI  ♦ DELTI 

246 

IF  (VI  .GE.  0.)  GO  TO  124 

247 

DELTI  ■ DELTI  - VI 

248 

CHNGI  ■ CHNGI  - VI 

249 

VI  ■ 0. 

250 

C*»  RE-EVALUATE  AT  NEM  POINT. 

251 

124  CALL  SETUP 

252 

CALL  DETNT  (N,  DET,  DETI) 

353 

TEMNR  ■ DELTA 

254 

TEMNI  > OELTI 

255 

TEMDR-VEL-DET 

256 

TEMDI-VELI-OETl 

257 

T EMDEN. T EMDR • T EMDR+T  EMDI *T EMD 1 

258 

IF  (TEMOEN  .EQ.  0.)  GO  TO  123 

259 

TEMRNU«TEMNR«TEMDR*TEMNI*TEMDI 

200 

ftM 1 NU>  TEMNI • T EMDR- TEMNR* TEMO I 

261 

C«*  EVALUATE  SLOPE  (RECIPROCAL  ACTUALLY  USED) 

• 

262 

SLR  ■ TEMRNU  / TEMOEN 

263 

SLI  ■ TEMINU  / TEMOEN 

264 

OELTA  * OET  * SLR  - OETI  * SLI 

265 

DELTI  - DET  • SLI  ♦ DETI  * SLR 

266 

IF  (K1  .EQ.  1)  PRINT  1330,  V,  VI,  DET, 

DETI, 

DELTA,  OELTI 

267 

C«*  CORRECT  PHASE  VELOCITY  TO  6EST  VALUE. 

268 

V ■ V ♦ OELTA 

269 

VI  ■ VI  «■  DELTI 

270 

TEMP  ■ V<*2  / ( TEMNR**2  ♦ TEMNI**2) 

271 

C**  MAS  INCREMENT  LARGE  ENOUGH  TO  PERMIT  EVALUATION  OF  SLOPE. 

272 

IF  (TEMP  .LT.  RLIM)  GO  TO  123 

273 

IF  (TEMP  .LT.  1.E34)  GO  TO  141 

274 

SLR  ■ SLOP 

275 

C*«  IF  NOT,  USE  EXTRAPOLATED  SLOPE. 

376 

SLI  - SLOPI 

277 

GO  TO  141 

276 

123  CONTINUE 

279 

C**  IF  SO,  FIND  1 - RATIO  OF  SLOPES. 

360 

▼EMDEN  * f«lR**2  ♦ SLI **2 ) 

281 

TEMDR  ■ SLR  • SLOP  ♦ SLI  • SLOPI  - TEMDEN 

262 

TEMDI  ■ SLR  • SLOPI  - SLI  • SLOP 

283 

TEMP  ■ (TEMDR**2  ♦ TEMDI* *2)  / TEMDEN**2 

264 

IF  (TEMP  .GT.  TLIM ) GO  TO  116 

2BS 

C*«  SLOPE  RATIO  TOO  GOOD.  DOUBLE  STEP. 

286 

141 

DP  ■ DP  * RATIO 

287 

GO  TO  117 

288 

118 

IF  (TEMP  .LT.  BLIM)  GO  TO  117 

289 

PRINT  130,  PK,V, VI, DET.DETI, SLR, SLI, TEMP, DBLOSS, NUMBER 

290 

130 

FORMAT  (1X,E14.6,E16.9,E13.7,6E10.3,I5) 

291 

C«*  SLOPE  RATIO  TOO  POOR.  HALVE  STEP. 

292 

IF  (NUMBER  .LT.  7)  GO  TO  126 

293 

PK  - PK  - DP 

294 

DP  * OP  / RATIO 

295 

INCA  * DP 

296 

JX  - JX  ♦ 1 

297 

IF  (K2  .EQ.  1)  PRINT  118,  PK,  V,  VI,  DET,  DETI 

298 

C**  STOP  ON  5 SUCCESSIVE  FAILURES.  MODE  IS  LOST. 

299 

IF  (JX  .LT.  5)  GO  TO  107 

300 

PRINT  810,  N,  FREQ 

301 

810 

FORMAT  (14,  G12.5) 

302 

3 

DO  801  I - 1 ,N 

303 

PRINT  800,  C(I),  Z(I),  GAMMA(l),  DPK(I),  RHO(I),  G(I) 

304 

800 

FORMAT  (10G12.5) 

305 

801 

CONTINUE 

306 

GO  TO  4 

307 

126 

PRINT  127  , N.  TEMP 

308 

127 

FORMAT  (7H  NUMBER,  13,  22H  FAILED,  SLOPE  RATIO  .F10.6) 

309 

C«*  UPDATE  ALL  QUANTITIES  FOR  NEXT  STEP. 

310 

117 

INCC  > INCB 

311 

INCB  ■ INCA 

312 

INCA  ■ DP 

313 

PV(3)  • PV( 2 ) 

314 

PVI (3)  > PV I ( 2 ) 

315 

PV( 2 ) ■ PV(1) 

318 

PVI (2)  - PVI ( 1 > 

317 

PV(1)  ■ V 

318 

pvi ( t > . v; 

319 

JX  ■ 0 

320 

DENOM  • V * V ♦ VI  • VI 

321 

LAMBDI  > -OMEGA  * VI  / DENOM 

322 

DB  LOSS  - -8686.  * LAMBOI 

323 

SR ( 3 ) ■ SR(2) 

324 

$R<2)  ■ SR( 1 ) 

325 

SR ( 1 ) > SLR 

326 

SI (3)  - SI (2) 

327 

SI  (2)  ■ SI(1) 

328 

SI ( 1 ) « SLI 

329 

GV  . V**2  / (V  - FREQ  * ( V-PV( 2 ) ) / INCB) 

330 

PRINT  118,  PK. V, VI .DET.DETI, SLR, SLI, TEMP, DBLOSS, GV, NUMB 

331 

118 

FORMAT  <E15.7.E16.9,E13.7,6E10.3.F11 .5,15) 

332 

NUMBER  « NUMBER  + 1 

333 

C**  CHECK  TOTAL  NUMBER  OF  STEPS. 

334 

IF  (NUMBER  .GT.  KO ) GO  TO  3 

335 

GO  TO  107 

336 

999 

STOP 

337 

1250 

FORMAT  (I3.8H  LAYERS ,, F10 . 1 , 3H  HZ) 

338 

1260 

FORMAT ( 6E1 0.4 ) 

339 

1270 

FORMAT  (8H  ATTEN  « ,G1 0.5, 5HDB/KM) 

340 

1280 

FORMAT  (8F14.5) 

341 

1320 

FORMAT  (/.6E18.9) 

99 
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